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ABSTRACT 


The object of this project was to define the problem of the assimilation 
of remote sensing data into mathematical models of atmospheric pollutant 
species. An object of remote sensing of the atmosphere is to enable recon- 
struction of the concentration distribution of trace species over a region 
based on the data available from the instrument. The data assimilation prob- 
lem is posed in terms of the matching of spatially integrated species burden 
measurements to the predicted three-dimensional concentration fields from 
atmospheric diffusion models. General conditions have been derived for the 
reconstructabi 1 ity of atmospheric concentration distributions from data typi- 
cal of remote sensing applications, and a computational algorithm (filter) 
for the processing of remote sensing data has been developed. 
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RECONSTRUCTION OF ATMOSPHERIC POLLUTANT 
CONCENTRATIONS FROM REMOTE SENSING DATA - 
AN APPLICATION OF DISTRIBUTED PARAMETER OESERVER THEORY" 
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ABSTRACT 

The reconstruction of a concentration distribution from spatially-aver- 
aged and noise-corrupted data is a central problem in processing atmospheric 
remote sensing data. Distributed parameter observer theory is used to de- 
velop reconstructibil i ty conditions for distributed parameter systems having 
measurements typical of those in remote sensing. The relation of the recon- 
structibil i ty condition to the stability of the distributed parameter obser- 
ver is demonstrated. The theory is applied to a variety of remote sensing 
situations, and it is found that those in which concentrations are measured 
as a function of altitude satisfy the conditions of distributed state 
reconstructibil ity. 
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I. INTRODUCTION 

In the remote sensing of tropospheric species, a ground-, aircraft-, or 
satellite-based platform performs an instantaneous scan of a region of the 
atmosphere and measures the species burden within the field of view. With air- 
craft or satellite remote sensing the platform is in motion and the fiela of 
view is constantly changing. An object of remote sensing of the atmosphere 
is to enable reconstruction of the concentration distribution of trace species 
over an entire region based on the data available from the instrument. 

The reconstruction of a concentration distribution from spatially-aver- 
aged and possibly noise-corrupted data is a central problem in processing 
remote sensing data. In the absence of a mathematical model describing the 
spatial and temporal concentration distributions, the reconstruction can be 
carried out by standard data interpolation methods. However, when a mathema- 
tical model exists, the problem becomes one of matching the remote sensing data 
to the model solution in such a way that the incomplete data can be used in 
conjunction with the model to produce an estimate of the region-wide concen- 
tration distribution. This problem of the matching or assimilation of remote 
sensing data into mathematical models for atmospheric constituents is the 
subject of this paper. 

There exist a few recent stjdies that assess the capabilities of remote 
sensing for monitcring regional air pollution episodes. For example, Barnes 
et al. [1] conducted a comparative analysis of satellite visible channel ima- 
gery in ground-based aerosol measurements. For three cases, each of which 
represented a significant pollution episode based on low surface visibility 
and high sulfate levels, the results show that the extent and transport of 
the haze pattern can be monitored from satellite data. The study demonstra- 
ted the potential of the satellite to monitor both magnitude and aerial extent 
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of pollution episodes. In a related study, Lyons, et al. [2] reported on a 
demonstration project showing that currently available synchronous satellite 
data can detect the aerial extent of large scale hazy air masses associated 
with sulfate and ozone episodes. 

A study related to that of the present work was reported by Diamonte, et 
al. [3] in which they considered the comparison of remote and in situ data on 
pollutant concentrations from point sources. They considered typical remote 
sensing geometries to provide insight on estimation of plume properties from 
these measurements. In a study also related to the present. Kibbler and 
Suttles [4] considered the estimation of unknown parameters in a pollutant 
dispersion model by comparing model predictions with remotely sensed air quality 
data. A ground-based sensor provided relative pollutant concentration measure- 
ments as a function of space and time. The measured data were compared with 
the dispersion model output through a numerical estimation procedure to yield 
parameter estimates that best fit the data. 

The object of this paper is to define the problem of the assimilation of 
atmospheric remote sensing data into mathematical models of pollutant behavior. 
Since the atmosphere is a three-dimensional system, models of pollutant behavior 
are of the distributed parameter type [5]. Remote sensing data represent spa- 
tial averages of concentrations, so that the assimilation problem is, in es- 
sence, one of distributed parameter state estimation. 

First, the concept of distributed state reconstruct!' bi 1 ity is developed 
for the class of problems of interest. That is, the first question to be 
faced is - can the desired spatial -temporal concentration distribution infor- 
mation be recovered from the measurements in the absence of noise. The deri- 
vation of general conditions that allow one to answer this question is the 
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subject of Section II. In Section III a variety of cornnon remote sensing 
measurement configurations and atmospheric models are tested for reconstruct!' 
bility. We conclude in Section IV with general observations concerning the 
inherent potential of remote sensing data in analyzing regional air pollution 
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II. RECONSTRUCTIBIL I TY AND OBSERVERS FOR DISTRIBUTED PARAMETER SYSTEMS 

Atmospheric pollutant models consist of partial differential equations, 
linear in the case in which the species does not react chemically or in which 
it is produced or destroyed by a first-order reaction of the form A This 
case represents a wide class of important situations and is the one to which 
we direct our attention here. Nonlinear distributed models must be handled 
by linearization and therefore also fall within the present framework. 

Our interest in this section is to derive distributed parameter observers 
for systems described by linear partial differential equations with inhomo- 
geneous boundary conditions characteristic of atmospheric models. An observer 
is an algorithm that processes measurements of the state of a system to yield 
an estimate of the entire system state. An observer is most frequently 
employed when not all of the states of a system are accessible for measure- 
ment. In the present application, we will be generally interested in only a 
single state variable, the measurements of which have limited spatial resolu- 
tion. The observer is stable if its estimated state converges to the true 
state after a sufficiently long time. The concept of state reconstructibility 
is useful as a condition for the stability of the observer. Thus, if a meas- 
urement strategy satisfies the condition of state reconstructibility, then 
the corresponding observer is stable, and, the state (i.e. the concentrations) 
can, in principle, be estimated from the measurements. The condition that 
allows the reconstruction of the system state on the entire field is called 
distributed state reconstructibility. Associated with distributed state re- 
constructibility, the concept of uniform n-mode reconstructibility can be 


-5- 


developed. Both conditions, n-mode and distributed state reconstructibility, 
will be applied, in Section III, to typical remote sensing measurement 
configurations . 

There exists some previous work on observer theory for distributed param- 
eter systems [6-8]. Kitamura et al . [6] formally extended the lumped param- 
eter observer to the distributed parameter case. Gressang and lamont [7] 
developed a more complete theory of the distributed parameter observer, includ- 
ing reduced order observers. An application of distributed parameter observer 
theory has been presented by Kohne [9]. The most complete treatment of observer 
theory is that of Dolecki and Russell [8]. In the current work, distributed 
parameter observers are derived in a form appropriate for application to the 
class of systems representing atmospheric species behavior. In addition, a 
result of the present work is an explicit relation between distributed parameter 
reconstructibility and the stability of the observer. Observer stability is 
demonstrated using a technique of Hale [10] in which Lyapunov stability theory 
is extended to function spaces. 

We consider the linear distributed parameter system, 

= L x u ( x , t ) + B(x,t)f(x,t) (1) 

defined for t > 0, x - D. The domain D is a connected subset of a d-dimensional 
Euclidean space with boundary surface 3D. The d-dimensional spatial coordi- 
nate vector is denoted by x. The state u(x,t) is a scalar function and L x 
is a linear partial differential operator with respect to x. It is assumed 
that the operator L y is well-postd. The input ;"(x,t) is a known scalar 
function and B(x,t) is a known coefficient. 
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The boundary condition on (1) is 

B u(x,t) = h ( x , t ) x € 3D (2) 

A 

where 6 x i s a linear, spatial differential operator of suitable order over 
3D and h(x,t) is a known function. The initial condition is assumed to be 
unknown or incompletely known. 

We are interested in considering three types of measurements: 

Case 1: Spatial ly-Independent Integral Measurement s 

The measurement takes the form 

w(t) = / H(x,t)u(x,t)dx (3) 

D 


where H(x,t) is a spatial weighting function. 

Case 2: Spatially-Continuous Measurements 

w(x,t) = C(x,t)u(x,t) (4) 

where C(x,t) is a square-integrable function, i.e., C € L^. 

Case 3: Spatially-Discrete Measurements 

w^t) = H.(t)u(x.,t) i = 1,2,... i (5) 

where w..(t) denotes a measurement at the i_th measurement location x^ . By 
taking the limit to small volumes of integration in (3), we can represent 
a system such as (5) by choosiny H(x,t) = ( t ) 6 ( x-x^ ) , i = 1,2,..., 1. 
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For the moment let us restrict the problem to one spatial dimension, i.e., 
0 < x < 1 . Accordingly, boundary condition ( 2 ) can be expressed as 


B Q u(x,t) = h Q (t) 

BjUfx.t) = h j ( t ) 


x = 0 
x = 1 


( 6 ) 


Then the solution of ( 1 ) and (6) with initial condition u(x, 0 ) = Uq(x) can be 
expressed in the form* 


u(x,t) = / $ (r,0;x, 

0 

rtA 


t)u Q (r )dr 


+ || $> (r,T ;x, 

Vo 


t)B(r,i )/(r,T )drdr 


n .-i 

0 0 


;x,t)g(r,T)drd* 


where 


g(x,t) = 2h j ( t ) 6 ( x- 1 ) - 2h Q (t)6(x) . 

★ 

The adjoint Green's function <t> (x,t;y,i) is governed by 


( 7 ) 


( 8 ) 


' j : * L%*(x,t;y , T ) = 0 


( 9 ) 


with the terminal condition 

* (x,t;y,t) = 6(x-y) (10) 

*The explicit form of operators L x> 6 q, and 6^ are assumed as follows: 

L x (‘) = a 2 (x,t) ^ ^ + a i( x » t ) ^97^ + ag(x.t)(*) 

6x 

6 0 (*) = a 2 (0,t) ^ + e 0 (t)(,) 

6j(*) = 0^(1, t) + 6 j ( t ) ( - ) 


and boundary conditions 
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V =0 


( 11 ) 


★ * 

= 0 

★ * * 

The operators L x> Bq, and Bj are the adjoints of the operators L x , B Q , and Bj 
respectively. 

The extension of the adjoint Green's functions to higher spatial dimen- 
sions is straightforward. In higher dimensions, (9) and (10) remain the 
same with the general boundary conditions 


B x ( t (x,t;y,i) = 0 x € 3D (12) 

★ ★ 
where B is the adjoint of the operator B . In general, we note that 4> is 

X X 

related to the Green's function $ associated with the system (1) with 

★ 

homogeneous boundary conditions by the relationship $(x,t;y,T) = <t> (y,x;x,t) 

The adjoint Green's function for well-posed distributed parameter sys- 
tems can be constructed in a variety of ways. Expansion in spatial eigenfunc 
tions and construction of the adjoint Green's function from eigenvalues and 

eigenfunctions is a powerful method for linear systems. Let us assume that 
★ 

L x has an infinite series of discrete eigenvalues { } , i = 1,2, Using 

standard methods, the adjoint Green's function that satisfies (9 )- ( 12) is 
found to be [11] 

* V* - X n (t ' T) 

<t (x,t;y,x) = 2 ^> $_( x )<t> n (y)e (13) 

n=l 

where the eigenfunctions { 4*^ ) , i = 1,2,..., are the solution of the equation, 
L x <t! i s satisfying the boundary conditions (11) or (12). 
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II.l Reconstructibillty Conditions 

The objective of an observer is to reconstruct the system state when 
the measurements are incomplete. To be able to reconstruct the state the 
observer must be asymptotically stable. 

An identity or non-reduced observer for the system (1) with measure- 
ments (4) takes the form 


3u(x,t) 

3t 


L x u(x,t) + B(x,t)/*(x,t) 

+ G[w(x,t) - C(x,t)u(x,t) ] 


(14) 


where u(x,t) is the observer output and G is a suitably chosen integral opera- 
tor with the kernel G(x,y,t). 

Before presenting a derivation of the observer, we will establish the 
conditions under which the system (1) and (4) is reconstructible. We define 
the reconstructibility kernel function by 


Q(x,y,t) = 


= fj 4*(x,t;r,T 


)C 2 (r,-rH*(y,t;r,T)drdT 


(15) 


It will be shown later that the observer (14) is stable if Q(x,y,t) has a so- 
called generalized inverse, i.e., if there exists P(x,y,t) such that 


/- 

D 


P(x,r,t)Q(r,y,t)dr = 6 ( x-y ) 


(16) 


By formal differentiation of (15) with respect to time and use of the 
properties of the adjoint Green's function (9) - (12), it is found that 
Q(x,y,t) satisfies the following Lyapunov equation. 
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■ - L*Q(x,y,t) - Q(x.y.t)L* ♦ C 2 (x,t)« (x-y) (17) 

with the initial condition 

Q(x,y,0) = 0 (18) 

and boundary conditions 

B*Q = 0, QB* = 0 (19) 

* y 

★ ★ 

where L^Q = QL^. Although Q(x,y,t) is formally defined by (l!i), it is impor- 
tant to note that Q(x,y,t) may be computed from ( 17 ) - ( 19 ) witnout using the 
adjoint Green's function. 

By using the identity 

ap (x,y.t) = - f f P(x,r,t) P(s,y,t)dsdr (20) 

D "D 

P(x,y,t) can be shown to obey the following Riccati equation * 

= L x P(x,y,t) + P(x,y,t)L 

P(x,r,t) C 2 (r,t)P(r,y,t)dr (21) 

with boundary conditions 

6 P = 0, P6 * 0 (22) 

x y 

P(x,y,t) may be considered as the kernel of the integral operator P defined as 

*The impact of observation error on the design of an observer can be assessed 
from (21) by comparing P to that from the corresponding distributed param- 
eter filter. 
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P/(x) = / P(x,y,t)/(y)dy (23) 

T) 

for f € l^. 

A linear distributed parameter system (1) and (2) with measurement (4) 
is said to be distributed state reconstvuctible if and only if Q(x,y,t) de- 
fined by (15) has a bounded generalized inverse P(x,y,t) for t > 0. It may 
be shown that Q(x,y,t) has a bounded generalized inverse when Q(x,y,t) is 
bounded and positive-definite for t > 0 [11] . The system (1), (2), and (4) 
will be defined to be unifcmlu n- r x?de reconstruatible if there exists posi- 
tive constants M^, and o such that 



<f> n (x)Q r (x,y,t)<j> n (y) dxdy < M, 


(24) 


★ 

for all t > 0, where ^(x) is the eigenfunction of L x and the modified re- 
constructibility kernel Q“(x,y,t) is defined by 


Q C (x,y,t) 


■/{ 


4 (x,t;r,x)C 2 (r,T)4 (y,t;r,x) drdi 


(25) 


The system is distributed state reconstructive if (24) is satisfied for each 
of the eigenfunctions. The uniform n-mode reconstruct! bility test (24) is 
useful when P(x,y,t) cannot be found directly from Q(x,y,t). Since it is 
straightforward to extend the concept of distributed state reconstruct!' bi 1 ity 

to measurement Cases 1 and 3, detailed discussion is omitted here. 

★ 

Positive-definiteness of the kernel Q implies that 



f , (x)Q(x,y,t)/(y ) dxdy > 0 


for all t > 0 and f - L^. 
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II. 2 Minimum Variance Observers 
II. 2.1 Observer for Case 1 

For the system described by (1) and (3), we define the reconstructi- 
bility kernel function by 


Q(x,y,t) 


■U 


l 


(x,t;r,T )H (r,i)dr / H(s,T)'t (y,t;s,T)dsdi 


(26) 


where Q(x,y,t) obeys 


9Q(x,y,t) 

at 


- L*Q(x,y wj - Q(x,y ,t)L + H (x,t)H(y,t) 

A y 


(27) 


with initial and boundary conditions given by \ ±8) and (19). Assuming that the 
system is distributed state reconstructible, the existence of the general- 
ized inverse P(x,y,t) of Q(x,y,t), that satisfies 


= L P(x,y,t) + P(x,y ,t)L 


at 


- f P(x,r,t)H (r,t)dr J H(s ,t)P(s ,y ,t)ds 
D D 


(28) 


will establish the observer for the system (1), (2) and (3). 

Following Meditch [ 12 ] , we define the cost functional associated with 


the observer as 


-13- 


i HI 


| J [u(x,0) - Uq(x) ] A x [u(y ,0) - u 0 (y))dx 


+ \ f f (.(t, - J~ H(x,t)u(x,t)dx^ 


dt 


( 29 ) 


where is an arbitrary final time, Uq(x) is an initial estimate of u(x,0), 
and 


A„ (•) 


■[/ 


1-1 


P Q (x,y){ ■ 


}dy 


(30) 


Pg(x,y) is a bounded, symnetric, and positive-definite weighting func- 
tion. The observer is found by selecting u(x,t) so as to minimize (29) sub- 
ject to (1) and (2). By minimizing the augmented functional, 


J = J, 


tf 

■n 


x(x.t) 


9u(x,t) 

at 


L.u(x,t) - B(x,t)/(x,t) ! dxdt 


>] 


(31) 


the result is the Euler-Lagrange equation, 


ax(x,t) 

at 


- L x >(x,t) 


H(x,t) 



H(y,t)u(y,t)dy 


(32) 


with the transversal ity conditions, 

X(x,0) = A x [u(y ,0) - u Q (y)] 
>>(x,t f ) = 0 


(33) 


Equations (32) and (33) constitute a two-point boundary value problem 
that may be solved by the sweep method. We assume the following Riccati 

A 

transformation for u(x,t), 
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u(x,t) = j P(x,y,t)A(y,t)dy + p(x,t) 


where the kernel P(x,y,t) and p(x,t) have to be determined. 
Substitution of (34) into (1), (32) and (33) yields 


= L x p(x,t) + B(x,t)/(x,t) 


j P(x,r,t)H (r,t)dr |w(t) - j H(s,t)p(s,t)dsJ 


p(x, 0 ) = u Q (x) 


B x p(x,t) = h(x,t) 


X t 30 


= L x P(x,y,t) + P(x,y,t)L y 

-j P(x,r,t)H (r,t)dr J H(s,t)P(s,y,t)ds 


P(x,y,0) = P Q (x,y) (39) 

6 /(x,y,t) = 0 , P(x,y,t)e y =0 x,y € 3D (40) 

A 

Equations (33) and ( 34 ) imply that p(x,t f ) = u(x,t f ) is the state esti- 
mate at an arbitrary final time t^. It is important to note that (38) is 
identical to (28). Thus we may conclude that the symmetric, positive-definite 
kernel P(x,y,t) completely characterizes the minimum variance observer. 

Equation (35) can be rewritten as 


L u(x,t) + B(x,t)/(x,t) 

X 


K(x,t) w(t) - / H(y,t)u(y,t)dy 

L D J 


1 
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where a time-varying observer gain K(x,t) is defined by 

K(x,t) =J P(x,y,t) H (y,t)dy (42) 

D 

The structure of the observer is identical to that of the distributed param- 
eter filter [13 ] . 

A 

We introduce the reconstruction or observer error e(x,t) = u(x,t) - u(x,t). 
Then we obtain the following equation for e(x,t), 

= L x e ( x »t) ‘ K ( x *t) / H(y,t)e(y,t)dy (43) 

OL A Q 

with initial and boundary conditions, e(x,0) = u Q (x) - u(x,0), and 6 x e(x,t) = 0. 
If the initial state is known exactly and the observer is initialized such 

A 

that u(x,0) = u(x,0), then the observer will reconstruct the state exactly. 

It is not reasonable, however, to expect that the initial state will be known 
exactly. It is, therefore, important to insure that if errors are present in 
the initial conditions applied to the observer that the estimate will converge 
to the true value of the state, i.e., the reconstruction error e(x,t) must 
have the property 1 im ||e(x,t) || = 0, for all e(x,0). 

t-KX> 

Asymptotic stability of the observer can be demonstrated by using (16), 
(26). (27), and (43). We will consider a Lyapunov function defined by 


V(e,t) 



e(x,t)Q(x,y,t)e(y,t)dydx . 


(44) 


It is first necessary to note that Q(x,y,t) is positive-definite and bounded 
from below. Then the time derivative of the Lyapunov function is calculated 
using (27) and (43). The result is 
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d_ 

dt 


V(e,t) 



e(x,t) H(x,t)H(y,t)e(y,t)dydx 


(45) 


0 D 

which is a negative-semi definite quadratic form. This is sufficient tc show 
that (43) is stable in the sense of Lyapunov [10]. 


II. 2. 2 Observer for Case 2 

In a similar manner to that of Case 1, we can obtain the minimum vari- 
ance observer for Case 2, i.e., for spatial ly-continuous measurements (4). 

The observer dynamics are described by 

A 

= L x u(x,t) + B(x a t)/(x,t) 


+ J G(x,y,t)[w(y,t) - C(y,t)u(y,t)]dy (46) 

D 

with initial and boundary conditions 

u(x,t) = u Q (x) (47) 

6 u(x,t) = h(x,t) , x € 3D (48) 

A 

where the optimal gain kernel G(x,y,t) is defined by 

G(x,y,t) = P(x,y,t) C(y.t). (49) 


The Riccati equation for P(x,y,t) in (49) is identical to (21) with boundary 

A 

conditions given by ( 22 ). The reconstruction error e(x,t) = u(x,t) - u(x,t) 
satisfies 

3e(x,t) 

at 


= F e(x,t) 


(50) 
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where the integro-di fferential operator F is defined by 


F e(x,t) = L x e(x,t) 



G(x,y,t)C(y,t)e(y,t)dy 


( 51 ) 


We can demonstrate the stability of the observer by using the reconstructi- 
bility kernel Q(x,y,t) defined by (15) and the Lyapunov function (44). 

Under the reconstructs 1 i ty assumption, the derivative of the Lyapunov func- 
tion becomes 


dt 


V(e,t) 


j e (x,t) C 2 (x,t)e(x,t)dx 
D 


which is a negative-semidefinite quadratic form. 


(52) 


11.2.3 Obs erver for Cas e J3 

For the spatially-di screte measurements (5), i.e.. Case 3, the observer 
is given by the following system: 


where 


and 


= L x u(x,t) + B(x.t)/*(x,t) 

£ 

+ S G.(x,t)[w.(t) - H (t)u(x ,t)J 
i = l 1 


G i (x,t) = P(x,x.,t)H.(t) 


= L P(x,y,t) ♦ P(x,y,t)L 


at 


£ 

- H P(x,x.,t)H (t)H (t)P(x ,y,t) 
i=l nil 


(53) 


(54) 


( 55 ) 
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with initial and boundary conditions given by (47), (48), (39), and (40). 

Under the distributed state reconstructibility assumption, the stability of 
the observer can be demonstrated. 

I I. 3 Conments 

The relationship has been established between distributed state recon- 
structibility and the existence of an observer. Distributed state reconstruc- 
tibility is defined through the existence of the generalized inverse to the 
reconstructibility kernel. The kernel associated with the observer gain 
satisfies the same Riccati equation as does the generalized inverse of the 
reconstructibility kernel. 

III. REMOTE SENSING MEASUREMENTS AND ATMOSPHERIC MODELS 

In this section we will test both the n-mode and state reconstructibility 
of common remote sensing measurements with models of atmospheric pollutant 
behavior. By far the predominant mode of remote sensing is to measure the inte- 
grated quantity (burden) of material between the ground and some known altitude. 
Thus, both cases we consider here involve vertically integrated data. Various 
assumptions concerning the horizontal characteristics of the measurements will 
be tested. Three-dimensional models of pollutant behavior are generally based 
on the atmospheric diffusion equation [5] that describes the flow and diffu- 
sion of species. The object of this section is to ascertain if the customary 
remote sensing measurements allow one, in principle, to reconstruct the detailed 
concentration distribution. The distributed parameter reconstructibility con- 
dition derived in Section II will therefore be tested in each case. 
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III.l Measurements in a Layer with Horizontal Homogeneity 

The vertical concentration distribution of a pollutant in a layer with 
horizontal homogeneity can be described by the one-dimensional diffusion 
equation. 


subject to 


3u _ v 9 d u 


(56) 

K i ■ h o< t > 

z = 0 

(57) 

k Tz- 0 

z = 1 

(58) 


where h Q is a given flux at the ground (z=0) and K is the turbulent diffusion 
coefficient. 

The adjoint Green's function for the system (56 ) - (58 ) is 

00 

2 

0 (z,t;z',i) = 1 + 2 ^ cos(mTz)cos(mrz' )e^ n7T ^ (59) 

n=l 


State reconstructibil ity is then to be assessed by condition (24) using the 
modified reconstructibil ity kernel (25). 

We consider each of the measurement types '3), (4) and (5). The condition 
for uniform n-mode reconstructibil i ty is (24), which is written for <t> n = cos(nTTz), 
n = 0, 1,2,. . . , as 

f 1 ~ 

0 < Mi < J J cos(n-nz)Q a (z,z' ,t)cos(nTz' )dzdz' < <°° (60) 

0 0 

For each of the three types of measurement, the integral in (60) is: 
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Case 1: Spatial ly-Independent Integral Measureme 

/ 2(nr ) 2 K(t-T ) j f u , x , x. 

e ' ' Ml H(r,T)cos(nTTr)dr 

t-0 0 

Case 2: Spatially-Continuous Measurements 

2(mr) 2 K(t-r) [ l L x , J 2 

I e v / |C(r,T)cos(nTTr) 

t-0 0 

Case 3: Spatial ly-Di s crete Meas ur ements 

fe 2 (n"»Mt-T)j^ „. {■> ) CO s(nitz j ) 

t-o li = l • 

for n = 0, 1,2, .. . . 

From (61 )- (63) , we see that uniform n-mode reconstructibility is com- 
pletely dependent on the form of the measurement weighting functions, H(z,t), 
C(z,t) and H.(t) and on the eigenfunction, cos(n*z). The condition (60) 

A 

implies that J H(z,t)cos(n~z)dz f 0. We may note that this inequality is 
o 

essentially equivalent to the observability condition derived by McGlothin 
[141. Similarly, (63) implies that the system state is reconstructible by 
point sensors if the sensors are not located at the zeros of any of the eigen- 
functions . 

In the remote sensing problem, the measurement weighting functions are 
often taken as H(z,t) = 1 or C(z,t) = 1. When H(z,t) = 1, the condition (60) 
holds only for n = 0 implying that the spatially-independent integral 


nts 


dT 


( 61 ) 


drdi 


(62) 


di 


(63) 
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measurements do not allow reconstruction of the system state on entire 

h 

fields. This can be directly checked by computing the reconstructibility 
kernel Q(z,z',t). The system with integral measurements cannot be distribu- 
ted state reconstructive since the generalized inverse of Q(z,z',t) * t does 
not exist, since Q is not an explicit function of z and z'. 

When the measurements are spatially-continuous and C(z,t) = 1, the system 
is distributed state reconstructive. From the definitions of Q(x,y,t) and 
P(x,y,t) in (15) and (16), we have 


oo 

Q(z ,z ' ,t) = t + -s- ^ -y cos(nTTz)cos(mTz' )( i 
tt K n=l n^ L 


e 2(mi) Kt _ l 


} 


(64) 


and 


00 

P(z,z',t) = 4^ 2 K ^ n 2 cos(n7rz)cos(nTrz' )/ i 

n=l L 


e 2(nr)^Kt _ j 


}- 


(65) 


We may note that the integral equation (16) is satisfied when it is recognized 
that 


★ . 

A mode associated with the eigenfunction (t> Q = 1 (n = 0) can be reconstructive 
and the appropriate observer is 


3u 

St 





u(z' ,t) dz 1 

- 


Stability of the observer can be demonstrated by constructing tne Lyapunov 
function 


V(e,t) 



e(z,t)te(z' ,t) dzdz 


l 
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00 

6 ( z-z ' ) B 1 2 ^ cos(n7Tz)cos(mTz' ) (66) 

n*l 


P(z,z',t) is bounded and Dositive-definite, and for t > 0, the series (65) 
is uniformly convergent. 


III. 2 Measurements of a Steady State, Point Source Plume 

The concentration distribution in a plume from a continuously emitting, 
elevated point source can be described by 


3u 

3t 


3 2 u 

7 


♦ K 


3 2 u 

7? 


(67) 


where t is the time an element of fluid spends in the plume from emission, 
equal to downwind distance x divided by the wind speed. The source is of 
strength q located at t - 0, y - 1/2, z « z H (0 < z H < 1). The boundary condi- 
tions on (67) are 


u(0,y,z) = q6(y - 

1/2)6 (z-z H ) 

(68) 

o 

II 

y = 0,1 

(69) 

K — = h 

z = 0 

(70) 

2 32 0 



!*• 0 

z = 1 

(71) 


The adjoint Green's function for this system is 


(WTJ K (t-T ) 


( ® vrvn ; k. 

My.z.t; y'.z'.x) =< 1 + 2 V cos(n7ry)cos(n7ry ' )e y 

L n=l 




cos(nrz)cos(nrz' )e 


(nm ) 2 K (t-T ) 


(72) 
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Consider first a scanning measurement performed at a horizontal position 


y = y*. 


w(t) 



J(z)u(t,y.,,z)dz 


(73) 


where the scanning data w(t) are taken on a coordinate that moves along the 
t-axis. o(z) in (73) is the altitude-dependent weighting 'unction for the 
measurements. When J = 1, the reconstruct bility kernel function becomes 


0 (t,y,z;y' ,z' ) = t 


4 ® “ cos(mryJcos(nTny*) 


K H ^ ifel <n 2 +m 2 )t . 2 


f (n^+m 2 )ir K„t 

cos ( nny ) cos (nny ' )< e - 1 


(74) 


The system is not distributed state reconstructible since the generalized 
inverse of Q(t, y, z;y' ,z' ) does not exist. Therefore, we conclude that 
the scanning measurement (73) cannot, in principle, allow reconstruction of 
the system state. 

The same results can be obtained for the following measurement systems: 


w(t) 


.1 

u(t,y,z)dydz 


w(t,y) 

w(t,z) 



u(t,y,z)dz 


u(t,y,z)dy 


(75) 


(76) 


(77) 
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In these cas^s, the reconstructibility kernel function Q(t,y,z;y ' ,z' ) cannot 
be written explicitly in terms of all the spatial variables y,z, and y',z'. 
Thus, the generalized inverse P(t, y,z;y ' ,z' ) does not exist, which allows 
reconstruction of the system state on the whole domain. As a rule, if 
Q(t,y,z;y' ,z* ) is expressible as an explicit function of all the spatial vari- 
ables and if it satisfies t‘ e uniform n-mode reconstructibility test, then the 
system is distributed state reconstructible. 

Indeed, we can show that the system state is distributed state recon- 
structible for the measurement 


w(t,y,z) = u(t,y,z) 


(78) 


In this case, we have 


i V I f 2(mrrK H t 

Q (t,y,z;y‘ ,z' ) = t + ir Z ? cos(niry)cos(mry' )< e -If 

n=l ( n-rr ) ^ L 

1 V 1 f2(mr)*Kyt 

+ -j- ? ~ cos(nrz )cos(mTz' Me -If 

K V n=l (nn) Z ^ J 

2{(mr) 2 K H +(nrn) 2 K v )t J 


♦*S I cos(nTTy)cos(mrz)cos(nTTy l )cos(nrTz' ) 

n=l m=l (n7T) 2 K H +(rmr)^K v 


( 79 ) 


The generalized inverse of (79) is given by 


” r 2(nn) 2 K H t VI 

P(t,y,z;/,z') = 4 K h 2. ("’> cos(m.y)cos(rmy' W e - 1] 

L n=l v 

v r 2 Cmr) 2 K v t Vl 

+ 4K W X (mr) 2 cos(n7Tz)cos(mTz’ ) \ e - M 

V n=*l ^ 

1 8 5 + (nm) 2 K v ^ cos(mTy)cos(imTz)cos(mTy' )cos(nmz ) 


2{ (n-rr) 2 K + 


where (80) satisfies the Riccati equation associated with the measurement (78) 
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IV. CONCLUSION 

This paper has examined the possibility of estimating atmospheric species 
concentration distributions from remote sensing data. Atmospheric concentra- 
tions can be modeled by partial differential equations of the diffusion type. 
Remote sensing data generally represent spatial averages of the concentrations, 
frequently in the vertical direction. The essential problem, therefore, is to 
assess the possibility of estimating the state of a distributed parameter sys- 
tem on the basis of spatially-averaged measurements. The theoretical basis 
of the assessment is a condition for state reconstructibil ity of distributed 
parameter systems. (The connection between state reconstructibil ity and the 
stability of the distributed parameter observer has also been developed.) 

A variety of remote sensing measurement configurations were tested for 
reconstructibil ity. It was found, not unexpectedly, that those measurements 
based on integration of tne vertical concentration distribution over the 
entire layer do not lead to distributed state reconstructibility, i.e., tnere 
does noc exist a generalized inverse of the reconstructibility matrix kernel 
and therefore do not afford the possibility of estimating the concentration 
distribution over entire field. Those measurement configurations that, on 
the other hand, enable sampling of the concentration at vertical positions 
lead to distributed state reconstructibility. 
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cmt runr fiorihmd parvwirr or b\ a unified approach 

bavr<J on rto U lenct-Hopf rtworx TW V 'opt oquanon for rtw 

cCTmjtion problem^ i\ #rmt4 vnitp (be «. esn manor error 

tnifnon l %in( rhr bout equation, fhrrr tyjts of thr optimal aoo dung 
Mimiion are nanwi«. fi\«4 point fnafintmaL anO fixed -lag 

wnoothrrx Finally. tbr rrujiti obtained are applied to etfimanon of anno 
ipberu wlfur concmfranom in the Totuxhima prefeemrr of 

Japan 


I Introduction 

A NUMBER of important physical phenomena ma> be 
modeled at- discrete-time distributed parameter sys- 
tems When estimation problems are encountered in such 
svsiems the measurements are also frequent!) discrete in 
time A gjeat deal of work has been earned out on estima- 
tion problem 1 for comtnuous-time distributed parameter 
system 1 [1)-|4! Tzafestas |5]. [6] and Nagamine n a! | 7 ] 
hast deroed optimal estimators for disoreie-lime distrib- 
uted parameter systems Tzafestas employed a Basesian 
approach where Sagamine er a! considered onh the filter- 
ing problem based on the Wiener-Hopf theory Recentl). 
Bencala and Seinfeld (?] have derived the optimal filter for 
continuous-time distributed parameter svsiems with dis- 
creie-time observations b> the Wiener- Hopf approach 
The object of this paper is twofold First. we seek to 
dense optima! filtering and smoothing ilgonthms for dis- 
crete-time dismbuted parameter systems b> a unified 
Wiener- Hor' approach Fixed-point. fixed-interval. and 
fived lag smoothers are considered Second, we wish to 
apply the results to the estimation of atmosphenc sulfur 
dioxide concentrations in the Tokushima prefecture of 
Japan 

II Description or the Distriilted Parameter 
Svstem 

Let i) be r bounded open domain of an r-dimensional 
Euclidean space with smooth boundary 3 D The spatial 
coordinate vector will be denoted bv x = (a,.- • -, a,) 6 D 
Consider a bnear distributed parameter system described 
by 

u( * -t- 1. x ) = t M y(k. x) + G(k . x )»(A, x). x € D 

(I) 

where ulk * 1. xi is an n-dimer.siona! vector function of 
the system »fi. z) is a vecto'-valued Gaussian process. £, 
is a bnear spatial matru differcnual operator, and G(A. z) 
is a known matnz function 
The initial and boundary conditions are given by 

e(0. x ) = e 0 (x) (2) 

r,irU S(* * l.t). ie»D (3) 

r<( j = <I J-(i-«<t»ai )'»• <<) 

w here « is an extenor normal vector lo the boundarv iD at 
» point t £ iD and o(£ i is a function of class r : on 3 D 
sansfving 0 « at { | « 1 S( * * 1. i I denotes a source func- 
tion at ihe boundary and is assumed to be known 
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Assume that w„( \ l is a Gaussian random function the 
mean and covariance functions of which are given b> 

£[*<><* >]=° (5) 

f[u 0 (a)ui(»)j * P 0 (a.>) (6) 

where £[ ] and the prime symbol denote the expectation 
and transpose operators, respectively 
Let Ute observed data be taken at m points, a'. - '. 
x" 6 D = D J iP and let an mn-dimensional cciumn vec- 
tor u m ( k | be defined b\ 

wjkl = Col[u(A. *').•••. i»( A. a")] (7) 

Let the observations be related to the states by 

x(k) = mk)ujk) + i{k). (81 

where :(k ) is a /> -dimensional observations vector at .he m 
observation^ points, a 1 .- • •. a“ e D. k ) is a known 
pxnm matrix, and t(k) is a p-dimensiona! vector-valued 
wnne Gaussian process Assume that the white Gaussian 
process wtk.a) in (li and r(/t) in (8i are statistically 
independeni of each other and also independent of the 
initial condition u 0 < » 1 Their mean and covariance func- 
tions are given by 

£[«(*. a)] =0. £[r(*)] = 0 (9) 

£[ w ( k a )» '( s. t )] = Pfk . a. t )6,,. a. * € D 

(10) 

£[fUlt(j»] = «(*>«*,. ( 11 ) 

where i,, is the Kronecker delta function and Q(k. x. i ) 
and Ri k i are symmetric positive-semidefinite and positive- 
definne matrices respectivelv 

111 DESCRIPTION OF THl ESTIVUTIOV PROBLEMS 

The general problem considered here is to find an esti- 
mate lit r. x k i of the state iiimi at time r based on the 
measurement data denoting a family of ci o i from 
o = 0 up to the present time k Specifically for t > k we • 
have the prediction problem for r - k the filtering prob- ' 
lem and for r < k the smoothing problem As in the 
Kalman -Buev approach an estimate ut t. a k i of ufr. x I 
is sought through a linear operation on the past and 
present observauon values as follows 
a 

u( T. X 'k | = 3 ^( T - *■ 0 )-( 0 I (12) 

where £i r. x. e ) is an nxp matrix kernel function 
To differentiate between the prediction filtering, and 
smoothing problems, we replace (12) with different nota- 


tion for each problem 
1) Prediction (t > It ) 

a 

u(r.x^k)- 2 -d(T. x.o)r(o) (13) 

• *o 

2 1 Filtering (r = * ) 

a 

fi(*.a/*) = 1 F\r. a.e)f(e) (U) 

3) Smoothing ( t < k ) 

a 

i(t.jr/k) = 2 T - * • e )?(• ) (15) 

• «0 

The estimation error is denoted by u(r. x fk ). 

fi(r. x/k ( = i»(t. *) - 6(t. x/) ) (16) 

The estimate u< r. x 'k i that minimise' 

Mu\- f[ i(r a k ) : ] (PI 


is said io be optimal mhere ! I denotes the Euclidian 

norm 
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Theorem } (Wiener- Hopf iheor>i A necessary and 
sufficient condition for the estimate £t r. x 'k Mo be opti- 
mal is that the follow mg_\*iener-Hopf equation holds for 
o * 0. 1. - ■ *. k and .t e D. 
x 

J £(t. j.o )£[:(» ):’(o)] = £[u( r. a );’(«)] 

••0 

( 18 ) 

Furthermore. (18) is equivalent to 

f[fi(r. x/k ).-'(«)] =0 (19) 

for a = 0. I. . k and a € D 

Proof Let F j( r. a. o I be an n\p matrix function and 
let < bt a scalar-valued parameter The trace of the covari- 
ance of the estimate. 

t 

w,(t.x/L) = J (F(T.a.o)-r tF±( r. x. o ))r(o ) 

• *o 


is given by 

•/(£,) = £ 


j(t. x I - ii( r. x/k ) - t ^ F^t. x.o )r( c 


= £ u(t. » k I : - 2 e£ 


j(t.x k) ^ F±( r. x. o );(o ) | 



^ £j( t. x. o l:(e I 

• = C 


A necessary and sufficient condition for £(r. x/k) to be 
optima! is that 

3 J( 6 , ) 

8t J 

that is 


0. 


£1 £'( r. A k ) ^ £j( T A. O ).-(0 ) 


for any n\p matrix FA r. a oi l sing the relation between 
the trace and inner product yields 


6 ( r. a A l ^ £jl r. a o 1:1 e ) 

• *0 


= IT 


61 t. x/k I ^ :'(o a. o ) 


= tr[£[u(T,A k ):'(o)]/j(t. a. o )] = 0 


Setting t. a. k ) = £[tf(r. x/k ):'(«>)] ui the above equa- 
tion. it follows that (19) is a necessary condition for 
tiiT. x k) to be optima! Sufficiency o' (19) also follows 
from the above equation Q.E.D 

Corollan 1 (Orthogonal projection lemma) The follow- 
ing orthogonalitv condition holds. 

£[filr. x /k)6({. v/k )] = 0. x. y 6 D (201 
where { is an\ ume mstant. for example. { < k. [ = * or 

f >k. 
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Proof Multiplying each side o( (l®i bs As. ». a) and 
summing from r*Oloe‘t yields 


lit v > 


a 

4 1 £ ;'<a)A(». x.a) 
• •0 


* 0 


Substituting (12 1 into the above equation yields (20i 

Q.ED 

Then the following lemma can be prosed 
Lmtma 1 (Uniqueness of the optima' kernel) Let 
At. a. e I and At t . > o I » \(t. a. o I be optimal matrix 
kernel functions satisfying the \4iener-Hopf equation (IKl 
Then it foils'* s that 


Mi.i.rliC • = 0_1 . 4 and x € £>. (21) 

and the optimal matrix kernel function At. a. b I is unique 
/Vos.i From <IKi we have 


A 

2 A(T.a o)£[.*(e); (a)] 

m • i 

* £[u(r.a).- l«l] 

A 

= ^ ( A t » . <i I * \(t a o ))£[(.-( o l:(o)] 

• •o 

Thus 

i 

^ V(t. a. • )fj;(o !;'(•)]* 0 

• •l* 

Multiply mg cash side of the abose equation b\ A' (t. a. at 
and summing from a = 0 to o = 4 yields 

A a 

2 ^ '(».a b)i'[.*(o !.*'(«)]'( t. a. a) = 0 

••o • • i* 

On the other hand from (ki and (III vs e hast 

£[ :i o I; l o l] = Wiel£[t.,'rli<..inl]Wi<xl _ /llo l£„ 

Then it fs'lls’sxs that 


^ ^ 'I* x e IW i o l£j u.l e lu„( o l] H ( o I ' t r x . o I 

0 ' 1 • • ( 

A 

• ^ h(r \ b IWi o lA( c tW l e I ' ( t. i . b I = 0 

• ■ i 

Since both terms on the right side of the abose equation 
are positise-xemidefimte because of the positise-definuenes' 
of Ate i a necessary and sufficient condition foi the aboxe 
equation to hold ts Vi t. a e 1 * 0. o = 0 1. .4 and x € 

f> Thus the ptsvf of the lemma is complete 0 E D 

In order to lactlitate the densanon of the optimal esti 
matsst' *e rewrite 1 1 8 • tn terms of the following coiollarx 
Coro/fun ’ The W iener Hopf equation ( IK) is rewtitten 
for the prediction filtenng and smoothing problems as 
follow s 

1 1 Prediction ( t > 4 ) 


A 

^ .41 t. x. b )£[;( o I; ’(b)] = £[w(t. a ).•'(«)] 

• •c 

( 22 ) 


for a = 0 1 .4 and i € D 

2) Filtering ( t = 4 ) 

A 

^ £( 4 . a « l£[ :( b ): ( a l] = £ [u( 4 . a );'(•)] (23) 

• »x> 

for a * 0 I .4 and a € P 

3) Smoothing 1 1 < 4 I 

i 

£ Al t. 4 a e )£[r(e )?'(•)] = £[•»( t. a I a l] 

■•r * 

(241 


for o * 0. I .4 and x € P 

In what follows let us den Me the estimation error coxan 
ance matrix function bs Pit. a. x 4 i 

P( t a i 4 I » £[ A( t. a 4 li ( t i 4 )] (2M 
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Ill Derivation or the Optimal Predictor 

In this section we derive the optimal prediction estima- 
tor by using the Wiener- Hopf theorv in the previous 

section 

Theorem 2 The optimal predicuon estimator is given by 
u(Z 4 \.x/k) = £,i(k. x/k) (26) 

r,6<* + I.</*)»S(* +!.{). tea/) (27) 
Proof From (22| and ( 1 ) we have 

i 

J Aik * 1. a. o )£ [c(e)c'(«)] = t,£'[w(*. a)c’(o)] 

• »0 

since ■>-(/. a I is independent of c(o). o = 0. 1. . Z From 

the Wiener- Hopf equation (23) for the optimal filtering 
problem we have 

i 

V ( y<( A -*• I. a. e) - t,P(Z. x. o ))£[;( o ).*'(<»)] = 0 

a * 0 

Defining S\ k. x. o ) by 

Ml« 1 . x. a ) = A(k -e I. x. o ) — t,£( k . x. o ). 

it is clear that i4(Z » I. a. o) + A’(Z * 1. a. o) also satis- 
fies the Wiener- Hopf equation (22) From the uniqueness 
of Aik l.a.o) by Lemma 1 it follows thai S(k * 
1. a. o ) ■ 0. that is. 

Aik - l.a.o) = t A flk.x.k). (28! 

Thus, from (13) and (I4i we have 

t 

fi( Z - 1. x/k ) = t, ^ £t*. a. c)r(o ) = i,6lk. x/k). 

a-0 

Since the forms of T ( and SiZ l.{) are known, the 
predicted estimate wlZ - 1. ( Zi also satisfies the same 
boundan condition as (3). Z -*• 1. 1 'k ) = Si Z - 1. 1 ). 
( € iD Thus the proof of the theorem o complete Q E D 
Theorem 3 The optima! prediction error covariance ma- 
trix 'unction Pi Z - 1. a. i Z i is given by 

PH - rl.i.t Z ) = f.PiZ. a. i Z)t; - 0(Z. a. » ). 

(29) 

T,P(Z * l.( y/k 1 = 0. ieiD (30) 

where 

C?tZ.a. i ) = GIZ.a)£>(Z.a. .t )G(Z..t ) (31) 

Proof From ( 1 ). ( 16). and i26i it follows that 
6(Z - 1. a /» I = i t u(k. x/k) G( Z. a )w(Z. a ) 

(32) 

and from (3). (16). and (27). 

r,u(Z -!,{/*;= 0. te»D (33) 

Then we have from (31) P(Z * l.a. v/k ) = £|6(Z ■* 
1. a Z >£(Z ■+ 1. y/k)} - i t Pik. a. y/k ILj * Qik. a. >) 
and from (33) EIT^Ik * l,(/k)£(k + 1. y/k)] = 
r,PtZ l.{. y/k) - 0 Thui the proof of the theorem is 
complete QE.D 

IV Deriv ation or thi Optimae Filter 

Let us derm the optical filter by using the Wiener- Hopf 
theorem for the filtering problem From (23) it follows that 

F\k * l.a.Z -s l)£[f(* -r l).- (a)] 

i 

* 1 Flk* l.a.o )£[.-(<) ).•<«)] 

• •0 

= £[«fZ -x 1. a ); (•)] (34) 

for o = 0. 1 . . Z - I 
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From (I) and the independence of c{ * ' and «•< A - 1.x). 
it follows that £|u(A 1. x )i'(o» = t',£ l“(A-X )-'(*» 

Applving the Wiener- Hopf equation (23) to the right side 
of the above equation yields 

£[w(A * 1. jr)r'(o)] = £. 2 £(*• x. a )£[x(a ).-'(•)] 

• •0 

(35) 


Furthermore, from (8) and the whiteness of u(A * 1) we 
hate 

£■[.•<* ♦ !).*'(■)] = W(A * I )£[*„( A - !).-'(■)]. 

Lei us introduce £.1 ) and | )i:i as follows. 

M ] o 


f.l ] = 


M ) 


(36) 


and 


[]t.=(L.|]| (37) 

Then from ( 1 1 and (7| it follows that 

yjk- )) = t.u.(A)~ wJA) (38) 

»JA ) = Col[C(A. a'MA.x 1 ). •••. 

(JoiA. a“)h(A.a-)] (39) 

Then we haxe for a < A 1. £|c( k - 1 ):'(o )] = W( A 

I)i.£|u„(A trial] Applying the Wiener-Hopf equation 
(23i to the nghi side of the abose equation yields 

£[--(A - l).* (o)] 

A 

= HU - 1)£. 2 £„U.e )£[.-< a );'(*)] 

«■* C 1 

(40) 

where 

f £|A. a'.o I 
TjA.o) = 1 

| £(A . x . o ) | 

Substituting (35i and ( 40 1 into (34i yields 

4 

2 Vjf A . x. a )£[:(<> )r'( a )] = 0. o * 0. A 

• «o 


(41) 


where 

KJk.x.o 1 = Fit- I a. A - 1 )HI A - Dt.FJA.ol 

- l 4 £( A . a . a ) -*• F( A ■» I.a.o) 

Since it is clear that £l A. a. a ) ■» Ajl A. x. a ) also satisfies 
the W iener- Hopf equation (23). it follows from Lemma 1 
that \ A t A. a o ) s 0 Thus we have the following lemma 
Lrmmc ’ The optima! matru kernel function £i A. a. o ) 
of the filter is given b> 

£( A ■* 1 . a . o ) = C, £( A. a . o ) 

-£(A * I.a.A * 1 )H{ A - I )£.£_( A . a ). 
o = 0 .l, .A ( 42 ) 


Theorem 4 The optima) filtering estimate Hk.x/k ) is 
given b> 

w( A ■* |. x/k ■* 1 ) = £,£( A. a 'k ) 

•+ £( A * I. a A * I )»( A ♦ 1) 

(431 

e( A - 1) = /(A * I) - HU * l)t.fiJA/A) 

(441 

£(0. a '0) = 0 (45) 

r,£(A - 1 ( A - I) = SI A * M) (€d£> (46) 

where 

6„(A A ) = Col[w( A . a 1 /A | . i( A . x*VA )] (47 1 
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Proof l ling 1 1 4 1 and (42) yields 
u ( k - i. */k * i) = m * I.*.* * ■):(* + I) 

2 f(*. *-•)*(•) 

•»0 

-£(* l.*,4 + l)tf<4 * DC. 

A 

• l FJk.'U(o) 

••o 

Again from (14) we have 

6{k + \.x/k ♦ l)= C.6( A. Jt/A ) 

+ £(4 * l.a. 4 l)r(4 -e I) 

Since we have no information at (he initial umc. it n> 
suitable to assume an miual value of u( k 1. x/k * 1 ) as 
i(Cl a '0> = £|w 0 < a)] = 0 Furthermore, since we know the 
exact forms of T ( and S{k + 1.4). the boundarv value 
u(4 * 14** I) also satisfies the same boundarv condi- 
tion as u< k - 1.4) Thus, we have T ( w< k * I. 4/ k ■* I ) = 
S(4 1.4). and the proof the of the theorem is 

complete Q.E.D 

Note that r (A •* I ) defined by (44) is rewritten by using 
the prediction value of (26) as follows 

r(4 - 1) = .•<* f I) - F(k * 1 )£.(* f 1 /A ) (4ft) 


or 

r( k - 1) = H(k - 1 )6Jk * l/k ) - c( A * 1) (49) 
where 

ijk - I A) = Col[i<4 - 1. x'/k ).•••. 

4(4 ♦ 1. *-/*)] (50) 

and 

i.U - I * 1 = u„(4 - 1) - 4„(A - 1. k I (51) 

►l A " il is termed the innovation process |8) |9] 

In order to find the optimal matnv kernel function 
fill - I. a. A - li lor the filtering problem we introduce 
the following notation 

FJ t. X i I = [pi t. a. a 1 /k ).•••. p(r. x.x m /k )] 

(521 


and 


r. .*' 4 I = 


Plt.a 1 4) 

P(t. a' k ) 
ph.x'.x'/k).-- .p{r.x'.x~/k) 

p{r.x m .x'/k). -.p{r.x m .x m /k) I 


(53) 


Note from the definitions of pjr. x/k ) and p m Ji/k ) that 
pjf. x/k)- £[i(r. x/k )i„( r/4 )] (>4) 

and 

p m Jf 4 ) = f [i„( r/A )£„( r/k )] (55 ) 

Furthermore we define the covanance mama of the in- 
novation process Kk — I) by T(* •+ 1/4). 

rU*l/*|«£[»U-HK(t+l)] (56) 

Then from (49 1 n follows thai 

T( A - 5/4) = H)k - I) p m Jk - 1/4) 

W (* - II 4 A(4 + I) (57) 


Then the following theorem holds 
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Thrartm The optimal filienng fun miitiv function 
hi * I a. * - It it pixen b\ 

F (A •» I . a . i * I ) * />_( A * I . ji >i ) 

«(* - nr-'u - i/aiosi 


ot 


m - i ji.* •* \)»fji - i.*/a> 

i(A * !/*)«(* ■* u/r'u + I) 
(59) 

where 

6(* + l/«)> (/* X(A 4- I) p m Ji 4 l/A))"(60) 
and 

XlA - J) * H\k * 1)A‘'(A * l)«(A * I) (61) 

Proof From the Wiener- Hopf equation (23) it fol- 
lows that 


FH - l a * * I )£[:(* 4 l)c(A + I)] 

A 

4 J Fy A - I. a. *)£[.-<« !;•(* - D] 

• •0 

= f[u(A - 1 a);'(A - II] 

Subs'iiutmf 1 42 1 into the aho\e equation yields 


F[k - l.a.4 •» II F 


:<A - II - H{k i I)U. 


I F.U n 1:(» I : (A - I) 

• ■0 


= £ 


u( A - 1 . j ) - L , ^ F[ k . i . 0 ).-( o ) } .•*( k ♦ I ) 

• » 0 


Substituting 1 14) into the n^ht side of the above equation 
and using t>» and the onhopona):t\ condition of (20» 
Meld' 


u« l 


1 ti - i, J fU a . * lr( a r ; *( A * 1 ) 

• I J 

= f[ii(i * 1 i k I; (A - H] 

= £[ill A - I « k tuJA - 1 1] W t A - 1| 
= r.u - la k \H (4 - II 


Lung the orthofonalitx condition of (20i fixes 
£[»U - I): (A - ll] ~ H[k - llfla.IA - I A |C\ 
Kk - I )] W ( A •» I) * X(A * I) 

= «|A - !)/>..( A - I 'A I 

H(k - I) i P[k * 1) 

* F(A 4 | A) (62) 

Then we have 

£U - I. a. A * l)HA - 1 A I 

= />.(* 4 1. a A )H (A 4 |) 

(631 


Thus < 51* i is demed In order to ihow the equixalence 
between t>hi and (59j. we use the followwf matrn in 
veruo.i lemrru. 


PH ( HPH -*)“«£(/• H R 'HP)"H R '. 

(64) 

From tSf i and (64 1 we hase 

F( A •» l a. A - || = f JA * l.a A) 

*<A - I A )« (A - |)A '(A * I) 

Then ( 59 1 it demed and the proof of the theorem is 
complete Q t D 

The equation for the optimal filtennf error coxanancr 
mains function p{k * l.a i A » 1 1 now must be de 
nsed 
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Theorem f> The optimal filiermg error cosanance matrix 
(unction plk I. a. v ■'* 4 I) it given by 
plk ♦ I. a. y/k ♦ I) ■ f(* * I. a. >/*) 

-pj* 4 l.x/k)HU 4 I) 

•r-'(* 4 |/*)W(* 4 I) 

*„<*-•• !.>'/*) (65) 

or 

r t + 1. a. y/k 4 I) = p[k 4 l.a, y/k ) 

-pjk 4 1. */* )£■(* 4 I/*) 

•*(* 4 l)p_(*4 | .y/k) 

«*> 

where 

p(0. a. > /0) = p 0 (a. y) (61) 

and 

T,p(* - l.(. \H 4 |) * o. ledD (Ml 
Proof From ( 1 1 and (43) we have 
6<A - 1 .x/k - l) = £(* - I .x/k) 

-FXk - l.a.* 4 1 1»( A 4 I) 

(Ml 

and from (3i and (46) 

r,£(* 4 l.(/*, 1)»0. (E9 D (701 

Using the independence propern between r(* * 1 1 and 
£( * •* 1. x/k I or fit * 4 1, i >k 1 yields from (69). 

p{k - l.a. t /* 4 |) = £[i(t * |.a * 4 |) 

■e ik * i. y/k - n] 

= />(*- I . a. i ,'A ) 

-/"(* - l.a.* - !>£[»■(* - I) 

r ( A - l)]f (* 4 |. , A - 1) 

-/■(A * l.a. A - 1 )HU - II 
£[uJA - I A l£( A - I. i A |] 

-£[£(* - 1. a A l£JA - I /k l] 

HU - I )n* 4 | ,y.k * I) 

Using (55 1 and (63) it follows that 

plk - 1. a. )/* 4 |) * p(A * |.a. y/k ) 

-pjk - I. x/k )H (A - I) 

Flk - l. i.A 4 |) 

= p(A 4 I. a. y/k ) 

-pjk 4 l.a A )H (A •> lfr~j ) 
fA * )/*)«, A 4 | )p m ( A ”, >/* | 

Thus (65) is demed The equivalence between (65) and 
(66) is easily shown by using (64) Since the initial value 
6(0. a/0) of 6U - 1. x/k ♦ 1 1 it lero from (45). it it clear 
thai /(0. a. s/0) = £|6(0. » /0)) = /pi a. >•) Multiplying / i f 

each tide of (70) by 6'(* 4 1, y/k 4 I) and taking the 
exgsectation yields T ( ?lk * 1. *. y/k 4 1) = 0. ( € 9D “ 1 f” 

Thus the proof of the theorem is complete Q E D 

Corollan J ijk 4 l/i* I) and pjk * I .x/k ■* l) 
satisfy the following relations. 

6. (A - I/* 4 |) 

= ij A 4 |/*) 4 FJk 4 |.* 4 |)r(A * |) 

(71) 

rjk * i.a 4 i) 

= p„.(* * I /*)£<* 4 |/*)H(* 4 ||**'(A 4 I) 

(72) 
or 

rj k 4 1 . a 4 1 ) > Fm J A 4 | /A 4 1 ) 

HU * l)*-'(A • I) (73) 

P.Jk - I, A - II = p m J A 4 l/k I — p m J A * | /A | 

■*(*4|/i)<|i4|);J*4|/l) (74| 
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or 

/>...< * - I *•»■> = p m J A - I /A IvU ♦ 1/A )• 

(751 

Proof From the definition* (411 and (JO) of Fj A •» 
I. A » 1) and 6„(A ■* 1/* ). it it clear that (71 1 . (7J>. and 
(74) hold From (60) and (74) it follows that 

/’„„<* - 1/A * 1) -F.J* - 1 '*:)♦(* * 1/*) 

{."(A* 1/*)- *<A + I) 

F..(*+ 1/*)} 

= F.JA ♦ l/*)(M* ♦ I/*) 

(/* *<A * I) p m Jk f I/*) 

-*l* - 1)F..(* ♦ 1/*)) 

= p.JA ■» 1 /A )d>(A l/* ) 

Thus. (7$) is denved and (73) is clear from (72) and (75) 

QED 

The present result corresponds to that of Santis er ol 
| • 7) which is an abstract form of the filter 
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V Drats ation or thi Eqi ations ro* Tm 
Optima: Smoothing Estimator 

In this section we dense the basic equa'ions for the 
optimal smoothing estimator bv using the Wiener- Hopf 
theory 

Lrnimc ? The optimal matrix kernel function B( r. k * 
I. a. o i of the smoothing estimator is gisen bv 

Jlr.l • 1. a . a ) * Bl r. A. a. e ) 

-£tvA - I . a. A ■» I )W( A - llt.FJA a). 
c = 0.1. .A (76) 

/’rocii From the Wiener-Hopl equation (24 » »e have 

< - i 

J Bl r , A - I . a o l£[ .•( e ).- ( a )] = £[ i/( r. a ): '( o )]. 


^ B\ t. A a . o |£[r( e ):'( a )] = £ [w( t. a lc'(o )]. 

• •o 

o = 0 . ■ • . A 

Subtracting the latter equation from the former yields 
file. A - 1 a. A - l|f[.*(A - l).-'(o)] 

i 

■* J ( El ». A * I. a. 8 I 

• «o 

-i(t. A. a. e ))£[.*( e ).-'(•)] * 0 
From (8i and (2?i we have 

£[.-(A - lie tal] = WlA - 1 )l.£[u„( A )c'( o )] 

A 

= tf<A - lit. V FJA.e) 

• •G 

•£[.*( e ): (8)] 

Then it follows mat 
< 

£ .V(t. A. a. e )£[c(e ).*'(o|] = 0 

• »r 

where 

fi'(T. A . a. e ) ■ itt.A-* I. a. e ) — E(r. A . a. a ) 

-Jlr.A * I. a A ■* ll«(A - llt.FJA 8) 

Since it is easilv seen that Bi r. A. a. e I ♦ Ni t. A a e| 
also satisfies the \A tenet- Hopf equation <24 1 from Lemrr.j 

I we have Mr A. a e ) m C and the proof ol the lemma is 

complete QED 

Theorr”- ' The optimal smoothing estimate ut r. a . A ■* 

I I is given b> 

i|t a A-I) = 6 |t a'A|-#|t.A-*1.a.A-I) 

►( A * II C7| 

r,il » ( 'A ♦ I ( * £( r ( I ( € iD (78 ( 


A = tt » | 



1 




Proof Frorr (IS) n follows that 
file. a * * I) * BIt.A * I. a. Jr * l);(A t- I) 
a 

* 2 B( r. * 1. a. e);(o ). 

••o 

Substituting (76) into the above equation vields 
fi( r. a/A t I) 

■ B(t. A +l,*,» 4 || 

jx(A t l>- W(A -e l)t. £ FJft.iM*)) 

i 

■* 2 ®< T • a.e)r(e). 

•*o 

Substituting (I4i and (1$) into the above equation yields 
£(r. a/A 1) * fi(r. a/A) 

* B( v. A -* I . a. A •» ) )r( A * I) 

Since we have no additional information about the 
boundarv value of u< t. a ). except for 5( r. ( ) and the exact 
form r ( . we have r { w< t. 4 'k * 1 1 = S(r. ( 1. 1 € 3D. and 
the proof of the theorem is complete Q E D 

Thrortm # The optimal smoothing gain matrix function 
B( t. A ■» 1. a. A * I) is given bv 

B(t.A - 1, a. A - I) * IJt. */k)ijg) 

17* 1)I"'(A * I /A) (79) 

or 

A * l.a. A ♦ 1) * J[r. x/k -4 I) 

H(k - l)*-'(A ♦ )) (801 

where 

J{ t x A » 1 ) = £„( t. a 'A )tl . 

(/ - AlA - 1 )/• f A * 1 'AH'* 


£.1 r » A I * [if t. x a 1 A ). . I(t. a. x" A )] 

( 8-1 


£(t.x.s A ) * £[fi( t. a 'A )i'(A. y fk )] (83) 

Pruvi From the Wiener- Hopf equation 624) it fol- 
lows that 

B(t A - I. x A * 1 If |r( A - 1).-(A - I)] 

A 

- V Bit. A - I. a. e )£[;(* ).•'( A * 1)] 

• ■ 0 

= f[u( t. a );'( A I)] 

Substituting (76) into the above equation yields 
Bit A - l.a. A - I )£[r( A -4 l)i(A f 1)] 

= £[ii( t. a/A )c'(* ■» D] 

On the other hand, from (48) and (49) we have 
£[r|A * l); (A - l|] = f[»(A - l)(»(A -4 I) 

■+ H(k * l)i.(A - 1/A))] 

= £(»(A -x 1).(A -r ))] 
-HA* I/A). (M) 

From ( 8 1 ami the independence of c(A ■» I) and 
fit t. x A | we have 

£[ilt. a/A ).*'(A ■* I)] 

* £ [ w ( t. a/A |fi a ( A * 1/A (A * 1) 

But from (26) and i)li it follows .nai 

6JA ♦ I/A) « t.iJA/*)- «I.<A) ( 8 ?) 

Then we have 

Bit. A - l.a. A * l|T(A •* 1/A) 

* i.|r.i/A):.H(A * I) 
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Let uv now derive the equation foi L(r. x. y/k 4 I). 
'■.'sing the orihogonalns condition (20i vields 

t(». a. y/k * I) * iC[ w« t . a)ii(A ■» I ,y/k * I)] 
'suhvtituling (69 1 into thr above equation vields 
it a. y/k •+ I) ■ 1(t. x. y/k )fj 

. Jr. x/ktL.HU * l)n* 4 I. }.k 4 n (16) 

F.\*i ■ and (78) it follows that ( A 4 ||t 0.(E 

c ■ Multi »In ing each tide bv ti'(* 4 1, y/k * |) and tak- 
ing i . v .^eilat.or yields 

r,E(T.|. i/|t||>0. { € »D (871 

Then the following theotem holds 

Theorem o Jt t » 'k * I) in (80) tt gisen bv 

J(r.x k - l)*J(f.*_OtWU- 1 k ) (881 

J(r x t) * />Jt, X/r) (891 

TfJir.l k - II = 0. t€ 3D (901 

/’fix.' From (8M and (59) n follows that 
LJ t > ( 11= LJ t. a )l. 

{/- A U * IW I* 4 | k) Fm Jk 4 | ()} 

Bui we have 

/ - Xt/ - rt\ 'r = *(/ - FR )'' 

((/ * PR)R ' - P) 

* (</ - PR\R-')'R 1 

= (*(*'' 4 F|)'' 

* (/ * jOm*'. 

Thus 

LJ ■ v k - II = f Jt « I.. 

(/ - R{k - I ) Fm J k * I A )) ' 
Therefore from i81 1 it follows thj'. 

Ar a * - II » L.lt.a k - I) (911 

and from (61 1 we base 
•fit. a k * I) = J(t. a k ).. 

(/ * R\k - l)p„(k 4 | /A))"' 

Then it follows that 

J[r x t|=(Jt » r) = pjr.n r| 

Since (90i is clear from (8 1 l and (9|). the proof of the 
theorem ts complete Q E D 

Let us now dense the equation for the optimal smooth 
mg error covanancr mams function pit. a. y k i defined 
b\ 

/>( t. a. t /k ) = E [d(r. x/k > * )] (92) 

From (7?) and (78i it follows that 
6 (t. a k ■* I) * 6(t. a k | 

-B(r.k 4 l.a.A + 1 >»( A 4 I) 

(9j) 

r ( i(t. ( 'k * I) = 0 c € »D (94) 

Then the following theorem holds 

Theorem ip The optimal smoothing error covanancr 
mams function />( t. a. > A 4 h ts given b> 

/>( r > 1 k * I ) =■ p( t, a. » k ) 

- LJ t. « k)C.H U 4 |^-j ) 

V + i >k )H\k 4 Dc.LJt » /k | 

(9M 

or 

/*(▼.*. i k ■* I) * p{*. 1 1 k ) 

- A r » k « |Q 

• *(k * I k | /»(* * 1 1 J t r . 1 ( II 

(96) 
(971 
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Proof' r.inii (W) it foUomr that 

p[T.x.y * - I) = p(r. x. y/k) 

+ Mir. It + l.a.A I )£[*{ A e \0 

*(A - 1 )] £ (t. A + I. t. A + 1) 

-B(r.k + J.jt. k * l)£[»(A + 1) 

■“'ir. >•/*)] - £[i(r. a/A )r'(A + 1)] 
■i'(r. k ■* 1.J.A + I). 

f[i(r. + I)] = £[ At r . x/k )i^( k /k )] 

* I) 

= ljt. t/k )i.HU * I) 

and 

£[r( k * l)i(r. »•/*)] = H(k ■» l)i.LJi.y/k) 
Thus we have 

p(r.x. y/k 1 ) = p ( r. x. y/k ) + B( t. k + 1. x. k + 1) 

m - i/a )*•(».*-*■ i. .».* i) 

— E( r . A + l.a.A + 1)W(A + I) 

•£.LJ r. i 'A ) - L m ( r. a/A ) 

•t.W(A -r l)£(r. A-t- l.y.k - 1). 
Substituting (79 1 into the above equation yields 
p{ t. i . v 'k - I) = pfr, x. yfk ) — £,( t, x 'k 

■Kk - IIT-'IA - 1 'A \H( A * 1 )£.£;( t.v/A). 
In order to denve (96). note that from (81). 

LJr x A = /( 7.x A - l)y ~'(A t- I/A) 
and from the matnx inversion lemma (64). 

H ' HFH - R\"H= (/- H R-'HPf'H JT'H 
Then we hase 


H (A - Ilf '(A - 1 A)«(A - I) 

= <<A - 1 A )/f(A * 1) 


and 

p(- . x t A - 11 = p(r. x. y/k 1 - /( r. x/k - I) 

•if'(A - 1/A )£( A -*■ 1)7(t. y/k* II 


MullipSing each side of (94) t> At r. y 'k -*• 1 ) and taking 
the expectation yields T £ />( t . y/k 1) = 0. { € iD 
Thus the proof of the theorem is complete Q.E D 

Corolmn 4 Jit. x k) sausfies the following relations 
J(r. x A - 11 = i4(r.a)/.(T + 1/A 1) (98) 


and 

A(r+ 1 . x /A ) = /)( r. a )J m { • /k ) (99) 

where 


/(r.a'/A ) 


JJ t 'k ) = 


( 100 ) 


[/(T.a“/*)J 

A(t. a) = fjr. a/r^ipjJ/r + 1/r) (101) 

A>(r. a) = pjr * 1. x/t)(p„(t/t)£;) (102) 


Proo f Letting ♦( A I ) be given by •( A + 1 ) = 
C.il £(A -» 1 )/•„„< A + 1/A))' 1 , from (88) and (89) it 
follows that J(t. x/k * 1 ) = p„( r. a/r )♦( r 1 >♦< t * 
2) ■ ■ ♦tA -*■ 1) and /„( r •* 1/A -f 1) = p mm ( r -f 1/r ■* 
I )♦( r 4 2) • • •( A + I ). From the above equations and 
(75) we have 

J(l. X/k » II = pj T. »/»)♦( T + I) 

p^lr * 1/T -*• ll/Jr * 1/A * I) 

= p„(r. a/T)t;i(>(T 1 /t)^''(t 1/r) 


P.llr-T \/t)JJt~ I/A 1) 

= >41 t. x )/.( t ■* 1 /A ♦ I ) 

From (88i and (89| it follows that /(r-» l x A) = 
p m ( r « I. a /» ♦ | )♦( r * 2 ^ ♦“< A ) and Jjr A I = 

F-.^' T ▼ >♦( ▼ — 1 )4>t ▼ — 2) *(A ) Thus, wi have from 

the abose eaualions 
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7(t* l.i A ) * p M { r - I . a r - I )( / - (It r - I ) 

p..(»- i '»»)( r— <» '/.(»/*) 

»)>■'.) 'j m i r/A) 

where the following equalnv derived from (66) has been 
used. 

p_( r* !.*/»♦ I ) * p^( r ♦ l.a/e) 

(/ /t(r ♦ \)p mm (* * l/t)) ( 10 J) 

Thus the proof of the coiollarv is complete Q E D 

Theorem II The optimal smoothing estimator is given 
b\ 


i(r. » A ) = £(r. a t ) - J 7(r.a//)#(/) 


/•*- 1 

(104) 

r ( fitr 4 A) * Jt».{). ( € »P 

(105) 

where 


HI) = tf(O*'(0r(/) 

(106) 


Furthermore the optimal smoothing error coxa nance ma 
ins function pi r. ». i A i is given bx 

ptr.a.i k ) * p(r. a. \ >i) - J h T • * 

fa . - I 7 —^ 

\l I - I )*(/)/ (T. y/l) (107) 
r,flr ( i t| = 0. ( € »D (108) 

Pro o' From I 1 '! and (80) 1 104 » can be directlx ob- 
tained and from t<»n (I0“) is dear Thus, the proof of the 
theorem is complete 0 E D 

\| St xtxnax O) THI Onixtxl SmisOTH|x.G 

Lsiixutors 

A F:\ed Part Smoother ( * - ft \tJ A = r — 1. r — 

: i 

Theorem !' The optimal fixed-point smoothing esttma 
lor is goer, bx 

■it ▼. a A -* 1 ) = £(*. a 1 I •* 7 (t.i A - 1 )8( A - I > 


(I0«i 

T .» k - l) = ytr.s A).. CIA -1*1 (110) 

ct* - i *) = (/-«(*- i)p„.u - uf ’ (iii) 

T. « T ) = p,( T. a/T ) (I ID 

r,£(r.t * - d = S(T.t). te »p (H3) 

r,JU a A - I) = 0. le»D (114) 


Furthermore the optima! fixed-point smoothing error co- 
vanance mains function pi r, a. > * - I) is gjxen bx 

p(t. a. * * - 1) = p(r. a. x A ) - 3(r. a A - l^*~^ 
r (A - I *)/}(* * l|T(t,v.-( - 11(113' 
r 4 P(r.(. y/k ♦ I) = 0. (CiP (116) 

B Ft ird Inte-xo! Smoothing Estimator (A = fixed, r = A 
-I.A-2 ) 

From Theorem 11 it follows that 
£(t-1.» A ) = £( t - I . a/t - 1 ) 

i 

- V y(r - 1 a 7)6(0 (117) 

fa.. ; 

and 

p( r - 1. a. i A)*p(T*l.a.i r — I ) 

t 

- J yn* I i /)£'•(/ /- I) 

/iif)/(T- i i o (ii8i 

Then the following theorem holds 
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Throrrm IS The optimal fixed-intentl smoothing ean 
maun is given h\ 

♦ I. t/k ) * 4(t ♦ I. a/» 4 I) 

•M(t * l.a)[fijt 4 2/4) - ij» 4 2/r 4 I)] (119) 
r,i(t* |.|/A)» J(,+ |.|). fP (120) 

Furthermore the optimal fiaedintersal smoothing error 
covariance matrix funeno: is given h> 

p(t 4 l.a, )/k) * p(r 4 |. a. >/» 4 |) 

->«(’ 4 1. a )(/>„„(’ * I/*) 
-p.Jr * 1 /t)M (t 4 1. i ) 

(121) 

r,p(» 4 i.|. »/*) - o. (eap (in 

Proof From (9Fi and (1I7| we have 
i|i4 1 . a 4 ) = i( t 4 1 , a ,'t 4 | ) 

4 

+ A(i - l.a) 2 Jjr - 1 '/)#(/) 

Bui from Theorem 11. 

A{t^2.x k ) * $(t «•* 2 . x '? ♦ 2) 

A 

4 I J( T 4 :,a /)»(/) 

!•>•) 

and from (4?i and i5^i 

ill ’4’. a ’4 21= it ’ - I.j ’4|) 

-ft’ 4 J.i.,4 2 lr( r * 21 

* m( • - 2 v ’ 4 | ) 

4 J{ • - 2. a ’ 4 21# t T 4 2) 

Thus »e hue 

i.t r » 2 * I - 6.1* ” 2 ’ 4)1 

A 

= 2 jjr 4 2 /)►(/» 

Then we have 

it *• l.a A ) = it ’ * 1 . a ’ * II ♦ d|t 4 |,a) 

[i.l’*2 A ) - i_( ’ * 2 't * ))J 

Frc>m (OSl and 1 1 1 K » 

pi ’ 4 ) a. i A I = /*( ’ * 1. a. t t * II 

4 

-*{’4 l.a) 2 JJf 4 2 

. 

*(/ I- \)R(DJJy * 2 /) 

/«!’ 4 I.t) 

From Theorem 1 1 

||t«2.i.i A ) = /*( ’ 4 2. a. > /t 4 2) 

4 

- J /(»♦). a//)t( * '(/// - l)A(/)J'(e * 2. y/l) 



*/ 
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y 


and from (M>| 

p ( t 2.x. jr/T f 2) = + 2. x. yf t ♦ I) 

t 2. x/t * 1 )<M » 2 /t + I ) 

■if(T 2)pJ» -f 2. >/t + 1). 

Taking into consideration that. from(103). J{r + 2. x/t + 
2) = />„(* * 2. x/t * 2/t •+ I) and 

;jT + 2/*)-;,.|f + 2/f + l) = 

- i JJ * + 2 //)* •'(///- i)J*(/)y;(rf 2/D. 

/■**2 

we have p{r * I. x. %/k) = p[r * 1.x. y/r ■* 1) - A(t 
* 2/A I - p m jr ♦ 2/r - 1. v) 

Since the boundary conditions ( 120 » and (122) are clear 
from ( 105) and ( 108). respectively the proof of the theorem 
is complete Q.E.D 

C Fixed las Smoothing Estimator (t = A+1.A = A + 1 
A. 1 = fixed I 

From Theorem 1 1 we have 
w ( k - 1. x,'k + I — A | = £( k ■* 1 . x 'k •+ | ) 

i-i-A 

* 2 J(k I . x .7 |r ( / ) ( 1231 
/**-: 

• plk 1 . x . v ‘k 1 A 1 
= p(A ♦ 1. x. i 'k ■» I) 



*[! I - HRi/l/fA - 1. t 7 1 

(1241 
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Then the following theorem holds 

Theorem 14 The optimal fixed-lag smoothing estimator 
is given h) 

i(t + l.«'» + lt A) 

«£,A(A.*/A-4A) + C(a.A*4 I.A) 

FJk 4 I + A/t ♦ | * A)r(4 t I -4- A) 

+ '*)£;)" 
•(i.|*/iiA|-4.(*/t)) (125) 

r,i(A - l.( A -4 I f A) * S(4 ♦ 1.(1. (£80 

(126) 

where 

C(a.4 * 1. A) = >4(4 + I. x )A m ( k + 1). 

• .-4J4 + A) (127) 

and 


/4J4) 


A[k .a') 
A(l.x m ) 


Furthermore, the optimal fixed-lag smoothing error covari- 
ance matrix function p(4 - l.a. i k * I •* A) is given 
b> 

p\k - l.a.t 4 * 1 - A) 

= p(4 * 1. a..»/A) - C(a. k * I. A) 

F m { k - I t A f * I 4 A)WU 4 1 4 A) 
r. J k * 1 - A k - A )C ( i.4* 1 . A ) - D(k. x) 
[/>„J4 4 ) - p mm t 4 4 * A )]/>'( 4 . »•) (121) 

r,p(4 * l.(. i '4 * 1 - A) = 0. ( £ 80 


(1291 


Proof From 1 43 • and (59| »e have 
x-14 - I. x 4 - l| = i,£(4.a 4 I 

■r7(4 ■» I. a/4 - l)»(4 * 1). 


From 1 12? i and the above equation it follows that 
£14 * l.a/4 -4 I * Al = £,£(4. a 'A I 


i-a 

* J 7(4 * l.a '4 )f (/) 

/•»- 1 

•4 7(4 -4 |.a// -4 1 -4 A)»(4 * I -4 A) 
From (88 1 it follows that 
7(4 * l.a 4 * I * A) = pjk -4 l.a/4 * I) 

•£W(4 * 2. 4 * I) 

CMk -4 3/4 -4 2) 

■£.V(4 + 2 -4 A/4 * I -4 A) 
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Substituting f-*5 > into the right tide yields 
/(A * l.*/* + I i- A) 

-pjk 4 \.x/k 4 l)£ji.-Ji(* 4 2/A 4 |) 

P m Jk 4 2/A 4 2)tMk 4 3/A 4 2) ••• 
tie (A 4 2 ♦ A/A + I * A) 

Repeating the tame procedure and uting (101) yields 
/(A 4 I, x/k 4 1 4 A) 

= A[k - l,’x)Ajk * 2) •• AJk 4 A) 

■p m J A * 1 * A/A 4 1 4 A). 

Thus we have 

7(A 4 l.a/A 4 1 4 A) 

A(A 4 1 4 A) = CU. A 4 I. A)p„JA 4 | 4 a/A * 1 4 A). 

(130) 

From |99i it follows lhai 
t*a t-a 

2 /(A 4 1, *//)#(/) = 2 />„(* + l.Jr/A) 

/* A — I /■*-! 

(*..(*/* ■»/ 

But from (29 1 we hase 

pjk * I. * 'k I = i,pj A. a 'A )L. 4 Q m ( A. a ). 

From (98 1 and (99) we have 
J(k.xrl) 

— A(k. x )/„( A - 1 '/) =pjT.*-T)t;(p„,jA/A)L;)' , /.(A,7). 

Then it follows that 
a - ^ i • ± 

2 /(A - 1. a /)!•(/) = £, 2 y(A. x/l)f[l) 

I* A - I /» A • I 

4 “ A - 1 

- A . a ) 2 [p m Jk A ) l .) JJ A / )A ( / 1 

/*i - 1 
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and 

*(A 4 i.x'A 4 I 4 A) 

* t,4(A.a/A * A) - Cfa.A 4 I. A) 
rjk 4 I * A/* 4 I 4 A >r(A * I * A) 

'i JJk/l)W) 

/•/• I 

Bui from Theorem 1 1 we h*\e 

ijk (4 AI-AJA'A)* 2 y„(A'/)f</> 

/»*- i 

Thu* we have (US) From (6$) and (124) ii follows that 
pik - 1. a. .* 'k I = p(k * I. a. y/k ) - J, - J : 
where 

a-a 

J, = I /(A - I. *//)<*'(///- l)rt(f) 

/•» - 1 

J(k - l..t /) 

y ; = JU 4 1.*/* •» 1 4 a )^ ” '( <r * I * a /A 4 a ) 
*(A 4 | 4 A) J(k * I. r/A - I * A) 

From (75i and (130) we have 

y. = C(x.k 4 I.AJp.JA 4 I 4 A/A 4 | ♦ A) 
ff(A * | 4 A )f mm ( A 4 1 4 A 'A 4 a ) 

C(i.A - I. A) 

= Cla.A - l.AIFJA 4 | 4 A A 4 | 4 A) 
HU-I* A)/vJA * | 4 A A * A) 

Cf.i.A - 1. A) 

Substituting (90 1 mtoy, nelds 

A - A 

y, = f>(A a I 2 yj a / »v’ *'(//- H 

- I 

/?(/)/.( A /ID (A i i 
But from Theorem 1 1 we have 

p( k . x i A ♦A|-f(li.J, i i A) 

» - a 

= - 3 y<A. a /K' '(/'/- l)Jt(/)y (A.>//> 

/■A - I 

and 

f.JA A - A) - p.JA 'A) 

A - ^ 

« - 2 yjA /k"(/ /- D*(/iyjA//). 

/*4- 1 

Then we have 

/>(A - |. x. y/k * 1 4 A) 

= />( A ♦ l.x. v./A)- C( a. A 4 |, a ) 

f,(A 4 | 4 A/A 4 I 4 A)W(A 4 | 4 A) 
p.jk 4 1 4 A/a 4 A )C'( > . A 4 1, A) - D(k.x) 
[PmJkk) - p m Jk/k 4 A )]/>'( A. j) 

Since tfct houndan conditions (126) and (120) are clear 
from (103i and t IOSi. respectively the proof of the theorem 
is complete Q E D 

Kells and Anderson 1 18] proved (hat the fixed-lag 
smoothing algorithm of Theorem 14 mas be unstable but 
Chirarattananon and Anderson |I9] denied a stable ver- 
sion of the algonthm It is possible to dense a compatible 
version here, although stabilus problems should not anse 
in our use of the algonthm of Theorem 14 as long as it is 
used over a finite time intern! 
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IX Amic mion to Estimation os Am 

PoLLlTlON 


Distributed parameter estimation theory hai recently 
been applied to simulated air polluuon data to demon- 
strate the capability of estimating atmospheric concentra- 
tion levels frcm routine monitoring data (10). (11) A 
problem identified in thesef early studies was how to 
specify the statistical properties of the assumed system and 
observation noise In this section we expand upon the prior 
studies in two respects First, we consider actual monitor- 
ing data for sulfur dioxide (SO] V in particular those mea- 
sured each hour during the period December 1-31. 1975 at 
four locations in Tokoshima Prefecture Japan (see Fig I) 
Second we apply the method of Sage and Husa [12) to 
estimate the unknown noise covariances in the system 
equation and measurements 

Hourly sulfur dioxide dau are available at the four 
locations shown in Fig 1 for the period December 1-31. 
1975 The data for day ft at location i may be denoted by 
e 4 (x'. r|. It is useful to average the data for December 
I - 30 to produce 

I * 

(eU'.r))*^ l r.(a'.r) (131) 

M A« I 


V 


Fig 1 


where we will consider December 31 as a day to test the 

algorithms 

If it can be assumed that the wind flows are such that 
there are no north-south sanations of concentration and 
that vertical mixing is rapid enough to eliminate sanations 
of concentration with altitude then the region can be 
considered to be one-dimensional alorg the east-west coor- 
dinate The SO ; concentration at any pamcular time can 
be assumed to be desenbed by the atmosphenc diffusion 
equation (13) 


f* l 

dr 


= a—, t Six. r ) 


(1321 


3 


where { is the wind selocity. a is a diffusion coefficient 
and S is the rate of emission of SO- as a function of 
location and time 

Equation bolds at any instant of time, but we 

desire an equation goseming the monthly mean concen- 
tration {() Although no such equation exists, we can 
formally aserage (!32i oser the 30 realisations (days) to 
produce 



One object will be to estimate the diffusion parameter a 
This parameter will in general vary with location and time 
of day although for simplicity we seek a constant value for 
the month Thus the first term on ;he right side of (133) 
becomes a ? : (u)/3r J We can form the residuals, u - c - 
(c) and : = r - (r). By subtracting (133) from (132) we 
obtain 




ic I 9c \ 9 : u 

‘8x " ( f 8a 


(134) 


Since wind dau are not available with which to evaluate 
the second and third terms on the left side of (134) let us 
rewrite (134) as 


iu 3 : w , . 


035) 


where «ia. / 1 includes those unknown features associated 
with the velocits terms 
The boundary conditions on (132) are 

1^ = 0. a = 0.1 (136) 

ei 

expressing the assumption that there is no diffusive flux of 
SO into or out of the region at the boundaries After 
averaging and forming the residual (136i becomes 

|^ = 0. a = 0.1 
dJ 


*/ 
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( 137 ) 



The problem is now to estimate irtx.il based on the 
data. 

x(x'.l) * »(x‘.») ♦ r,(r). !■ I.1.J.4, (138) 

Since hourl> data are available. (I3S) can be cast into the 
discrete-time form 1 1 ). 

y{k * I. x) *£,«(*. x) ♦ w(x.l) (139) 

with t,» I ♦ad’/Jjt 3 Observation error is estimated 
from the mean square error of predicted values and ob- 
served da la 

2 (:,(*> -!,(*/*- D). i*l. 2.3.4 

<■ i 

(140) 

An index of overall estimation error is 

a 

J=2FA{>) (141) 

• ■ I 
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To apply discrete-time distributed parameter estimation 
theon to predict air pollution levels, we must consider 
thiee problems The first problem is now to simulate the 
distributed parameter system The second is how to de- 
termine the covariances of system and observation noise 
The last is how to determine the diffusion coefficient a For 
the first problem we use the Fourier expansion method and 
approximate the original distributed parameter system by a 
finite-dimensional system For the second problem, we 
apply the algonthm of Sage and Husa [12] that necessitates 
the simultaneous application of the optimal filtering and 
smoothing algorithms For the third problem we apply the 
maximum likelihood approach in the smoothing form [14] 
We now consider these problems in more detail 

Founrr Expansion Method It is well-known that the 
state u< A. xl of the distributed parameter system (I39| 
with boundarv condition (13^1 can be represented by using 
the eigenfunctions e.t x ) as follows. 

w( A . x ) = 2 »,( A )♦,(*) (142) 

.» t 


where 


and 


£>(x) = A, ♦.<*). x £'(0. 1 ) 


aoui 

a* 


= o. 


{ = 0.1 


LJ 



(143) 


Au /Ed placement of eq no’’ 


( o ( « )d,(x )di = t,j 
J o 


A, is the eigenvalue of t, corresponding to ♦.(x) In this 
case, it is easily seen that the eigenfunction ♦,( x ) and the 
eigenvalue A, are given by 

♦ ,(*)=!. #,(x) = »2 cose.a. t = 


and 

A, = I - op 1 (i - 1) : . i = 1.2. ••• f"**) 

Then it*. x /k ). pi r. A /k ). and A(r. x) can be rep- 
resented as follows * 

m 

fi(r. x/*) = 2 i, y (»/* )♦,(*) 

I* I 

p(r. x. y/k ) = 2 A y ( )♦,(* )♦,’(> ) 

»./■ ' 

4( t. a ) = 2 «»< T )4,(x) (US) 
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Lei u* approximate the»e infinite expansion* b* the first A' 
term* and define the following matrices and vectors. 

i(r/0 - Col [-,<»/*). ••.*,<»/* )]. 

Mr) e Cot [«, I f). •••■**(»)], 

A = diag[A„- A,). 


and 


K t/A)' 


QU) 


♦ = 


#ari( »/*)■• •.#*»<*/*) 

♦l|(*). • * •.*!*(* ) 


*v, (A)-'. «**<*). 


where q, t 4 I denote.* the (i. j Ith Founer coefficient of 

(?U *. .») 

Then, from Theorems 3-5 we have 
i( 4 * l/k* l) = A i(4/*)f F(k f 1 )r( A - I) 

F(4 * 1) * f»(4 * 1/4 )♦ « (* * I) 

(«(4 - n»m * 1 /k)*fi) 
•f4* l)-r*(*fl)]'. 

P[ k - J/*) * A P(k 'k )A' + 0(4). 

P(k - l/k -*•))=(/- Hk - l)/f(4 + I)*) 

■P(k*\/k) (146) 


Furthermore from Theorem J2,we base 
lit t — 1 k ) = u( t - ) ■£- I ) ♦ A{ t I )♦( 6( r - 2 -4 1 

-6(t- : r - i)). 

<4( t f 1) = />|T- It - l)AF''(: - !/>)♦•'. 

P ( T •* | 4 I = Fl T I T t i) 

~A( r- I >♦(/*(▼ - )/4 ) 

-F(T*l/t|)*4'(t*l) (i<n) 

Note tnat the fixed-interval smoothing estimator does not 
depend on the matrix ♦ which reflects the effect of sensor 
location 

Petrnvination of ihr Sour Covariances In order to 
determine the unknown covanance matrices of the system 
and observation noises, we adopt Sage and Husa's algo- 
rithm 1 12) given by 

i 1 

0U> = 7 l (i(;/4l- Ai(;- 1/4)) 

y*i 

■ (6{j/k) - Au{ j - \/k)Y (148) 

and 

*<*> = { 1 ( 4 ( 7 ) -#(;)♦£(;/*)) 

i-i 

(*(/>- H(j)*6(j/k)Y (149) 

where Qt 4 1 and X < 4 ) denote?' the estimated values of 
Qi k 1 and At 4 j. respectively Note that in these idenufica- 
uon algorithm the fixed-interval smoothing estimate u{j/k) 
is used 

identification of the j inkntm-n Parameter a To determine 
the unknown parameter a we use the maximum likelihood 
approach in smoothing form 1 14J The log-likelihood func- 
tion r< 4 mis given from 1 14) by 

r(4.«) = ♦ Y* , I 
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* -Af In (2* I - 2 In do I ,( j;j - l a ) 


ORIGINAL PAGE IS 
OF POOR QUALITY 


TW, * - E (r'( /. a )#■'►( y. « ) *• [Hj'k.fl) 

/• I 

-Hj- I k.o)\Q-' 

(tily A. a) - Hj - I/A . « )) ) 
r[j. a I * :()) - Nijfyj/J - I. *1 { 

I.f / y - 1 a I ■ f[r( i. *)]. 

where p is the dimens.. of ;<Ai and in j k - l.ci 
denote- u i * - 1 1 under the condmor. that the unknown 
parameter is assumed to be u 
To maximire Yik.a) »e use the following gradient 
method 

a . t - a - Cl.' I> ( A. a, I 

?>( k a l 

a. I = — ; (1511 


>,U- 


iu 
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where Cut is a suitable matns Therefore we adopt the 
following recursive algorithm to idennfv the unknown 
pa'ameters Q. R and a 

1 1 Make an initial guev» a, of a 

2i Compute - and - 1 hv using ( 145 1 and tl4<3| 

?i Compute a t" using il*l'i 

4 compute (*' c i and Rid t h\ using 1 145 1 and 1 140 1 
5 1 Rerurrv,three hv changing / to i - I and repeal until £ 

these value- do nc»t change 

'u» ien<i>.' Ren. \k e uve the observed data from De- 

cembr I- 5, ideftv the unknown parameter c and 
noi-r covanan-e- C‘ and R Mie - fou' iteration' the alg.- 
ruhn. for determining a converged to the value <J = Ota. I 
The Fourier e\pan«ior ha> beer, truncated at ' = 4 The 
e't mated diagonal elements of noise covarun-c matrices 
arc 

Q R = o :« 

(>.. = 1 4v R ;; * on 

£>„ = 5 '5 R„ * 1 <* 

JM R„ - I 54 

To conside’ the effect o' the nurr£’ and location of A [ 

monitoring stations we assume that we have data at onlv 
one monitoring station In iho case from ihe previous 
results of Kumar and Seinfeld |l‘] and Omalu er al 1 16) 
we expect that the optimal sensor location is closest to the 
houndan Thus either a 1 or x* is Ihe optimal single sensor 
location among the four monitoring stations, a'. x : . a \ «* 

In Table I we shovs the values of R4til and J for several 
monitoring station' 5k e sec lhat Ai:umi or Matsushige is 
optima! for the one-point sensor location case Similar 
conclusion* hold for two or three moruionng stations 
Finalls we illustrate ihe actual observation data and one- 
hour ahead predicted values for December 31 in Figs 2-5 
for Ailumi kitaj^m* Kawauch. and Matsushige. respec- 
tive!) Table 1 Figs 2-5 

Compj’non with Oihr ■ App'itoche' li is o' interest to 
compare results of the present filtering and smoothing 
approaches with others available for air grollutiot esuma 
tior \ke consider therefo-e the same SO ; estimation 
rri'blem hv the following methods It ^R-model 2 > per 
sistrnce and 3i weighted ensemble 
The 4R-modei method is based on the following 4R< p i 
mode! 


t 


1 . 2.3 4 




(1521 


where the u\ s nre the concentration levels ai time * and 
at morvionng station x’.c l .a : . a t arc the cone 
spons 4 itig A ^-parameters and the c, s are re*idual' kke 
usesi the Levinson algorithm to determine the AH* 
parameters while the optimal order /> of the .4R-process is 
determined bv using the minimum Akaikcs information 
criterion (AlCi |20) Then the one-hour ahead predicted 
concentration is given b> 


w . i - 1 




(153) 


and the piediction error variance is 

• lit' 




■i l 


— ^ ( u"' - u"‘ 
14 on “* • 

* 1 . I 


0341 


Table 11 show* the .^-parameters and minimum AlC 
saluc at each monitoring station 

The persistence method consists mereK of using the 
observation data v', , as the one-hc>ur ahead prediction 
value u, , . , 

The weighted ensemble method uses the mean of the 
past observation data at each time k weighted bv a linear 
function of the source strength as the prediction value at 
lime k Based on the number of emission sources the 
weighting functions are assumed here to he 0 15. 0 41. 0 26 
and 0 lb at a’. a’, a 4 , respecltselv Table III shows the 

performance criteria of the four methods From Table 111 
we can see that the present method possesses almost the 
same accuracv as the AR-model method Bv multiplving 
each eigenfunction coefficient bv the corresponding eigen 
function and summing them however the present methsnJ 
enable- us u estimate concentration' over the entire re- 
gi r. There'orc the present method is more powerful than 
the CR-mode) method 
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Opt.ma; estimators for dis.retr tmr distributed parame- 
ter svstems ha'e beer densed based on \k tenet- Hc>pf 
theors A ns'ta'-’c point of the present work is that the 
smooth ng estimat rs hase beer dtrised bs the same ap 
proach a- the filter thus providing a unified approach for 
ihi- cIjs- of distributed parameter estimation problems 
The estimation algorithms have been applied to the prob- 
lem of predicting atmospheric sulfur diovide levels in the 
Tokushima prefecture of Japan 
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ESTIMATION OF ATMOSPHERIC SPECIES CONCENTRATIONS 
FROM REMOTE SENSING DATA* 


Sigeru Omatu and John H. Seinfeld 
Department of Chemical Engineering 
California Institute of Technology 
Pasadena, California 91125 


ABSTRACT 

A basic problem in the interpretation of atmospheric remote sensing data 
is to estimate species concentration distributions. Typical remote sensing 
data involve a field of view that moves across the region and represent inte- 
grated species burdens from the ground to the altitude of the instrument. 

The estimation problem arising fro" - this special measurement configuration is 
solved based on the partial differential equation for atmospheric diffusion 
and Wiener-Hopf theory. The estimation of the concentration distribution 
downwind of a hypothetical continuous, ground-level source of pollutants is 
studied numerically. 
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I. Introduction 

In the remote sensing of atmospheric species, a ground-, aircraft-, or 
satellite-based platform scans a region of the atmosphere and measures the 
species burden within the field of view. An object of atmospheric remote sen- 
sing is to reconstruct species concentration distributions over a region based 
on the data available from the instrument. 

There exist two recent studies that assess the capabilities of remote 
sensing for monitoring regional air pollution episdoes (1,2]. Diamonte et al . 

[3] developed theoretical results for the estimation of point source plume dis- 
persion parameters from remote sensinq ’’ta. In a similar vein. Kibbler and 
Suttles [A] studied the estimation of unknown parameters in a pollutant disper- 
sion model by comparing model predictions with remotely sensed data. No results 
have yet been reported in which actual remote sensing data have been used to 
estimate species concentration distributions. 

The present paper deals with the theoretical foundation of estimating atmo- 
spheric concentration distributions from remote sensing data. Since the atmosphere 
is a three-dimensional system, mathematical models of pollutant behavior are of the 
distributed parameter type [5]. Remote sensinq data usually represent spatial 
averages of concentrations, so that the estimation problem concerns a distribu- 
ted parameter system with spatially integrated, scanning data. Although dis- 
tributed parameter state estimation has been considered extensively (see, for 
example, [6] and [7]), such problems with scanning and spatially integrated 
measurements have not been considered previously. The purpose of the present 
paper is to derive the required optimal estimators for the scanning and spatially 
integrated measurement case by a unified method based on the Wiener-Hopf theory. 
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In Section II, we define the remote sensing data analysis problem mathema- 
tically. Sections 1 1 1 -V I are devoted to derivation of the optimal prediction, 
filtering and smoothing algorithms for the problem by Weiner-Hopf theory. Fi- 
nally, in Section VII we present a detailed numerical example of estimating the 
concentration distribution downwind of a continuous, ground-level line source to 
illustrate the application of the thoery. 

1 1 . Problem Statement 

We consider a single atmospheric species (nonreactive), the mean concen- 
tration u(t,Xj ,X 2 ,x.) of which over a certain region is described by the follow- 

★ 

ing form of the atmospheric diffusion equation [5], 


££ + V £i!_ + v — = — I K ( x ) 


oU 


3u _ 3 


3u 


3t 


1 9x 


2 3x. 


ax. 


V 3' ax 


3 


w(t,Xj .x^.x^) 


( 1 ) 


where Vj and Vr, are the mean velocities in the x^- and x^-directions , respec- 
tively, K v (x^) is the vertical turbulent eddy diffusivity, and w(t,Xj ,Xp,x^) is 
a randor disturbance accounting for inaccuracies inherent in the basic model. 
The initial condition for (1) is u(t .Xj.XgfX^) = u q (x 1 , and typical 

boundary conditions are 


M x 3 * 3 Xj = S(t,Xj ,* 2 ) , 


x 3 * 0 


-u 

ax. 


= o, 


x 3 = h 


( 2 ) 


where Slt.x^x^) is the ground-level species source emission rate, presumably 
a known function, and h denotes the upper vertical boundary of the pollutant- 
containing region, for example, the base of an inversion (stable) layer. For 
convenience, we denote the coordinate vector by x and let 


L 
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] = - V 


Mil 

1 3xj 


V Mil + _M 

2 5X2 9Xg 


K (x.) 
, v 3 


Mil\ 

ax 3 I 


In this form of the atmospheric diffusion equation, turbulent diffusion in the 
horizontal direction is neglected relative to transport by the mean flow, a 
common assumption in treating atmospheric diffusion problems [5]. 
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Assume that the concentration of a species over a fixed spatial domain D 
with its boundary 3D is of interest. Let us define the operator r^, C € 3D 
as follows. 


* r 


- K (x,) 
v 3 3x^ 

3[ • ] 

3x, ’ 


x 3 = 0 


x 3 = h. 


Let S(t,t) be 


S(t,0 - 


S( t,Xj pXg) » 
0 


x 3 = 0 
x 3 = h. 


Thus, (1) can be represented as 


■ L u(t.x) + w(t.x) (3) 

C*A A 

and (2) can be written as 

r r u(t,0 = S(t,£), C € 3D. (4) 

We assume that the initial condition u q (x) can be represented as a Gaussian 
process with statistics, 

E[Uq(x)] = u Q ( x ) 

E[(u 0 (x) - u 0 (x)J (u 0 (y) - u Q (y)) ] = P Q (x,y) (5) 

and the random disturbance w(t,x) is stochastically independent of u q (x) and is 
a white Gaussian process with statistics, 


Etw(t.x)] = 0 
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( 6 ) 


E[w(t,x)w(s,y)] = Q(t,x,y)5(t-s) . 

We assume that the remote sensing measurements are taken at time t^ over 
a view volume D(k) consisting of M pixels, as shown in Fig. 1. Since the sens- 
ing platform may be in motion, the field of view, in general, moves with time 
across the entire spatial domain D. We assume that the shape and extent of the 
field of view D(k) remain fixed and only the location of the centroid of each 
pixel changes with time. The ground-level location of the centroid of each 
pixel of D(k) is denoted as ( x^ 1 ' k ^ , x^^.O), m = 1,2,..., M 

We are interested in considering the vertically integrated measurement 
given by 


Z m(k) (t k ,n) = / 0 ^m(k) (x 3 )u(t k’ x l (k),x 2 ( ),X 3 ) dx 3 


+ v ( t x m < k ) x m ^ k ) h ' 

V 'V 1 ’ 2 ’ V 


(7) 


m = 1,2,..., M, n = 1,2,..., N, h < h 2 < ... < h N 
k = 1,2,... 

where ^ ( x^) is an a 1 1 i tude- dependent i nstrument weighting function, and h^ is the 
vertical position of the scanning sensor. Physically, Z m ^)(t k ,n) represents 
the vertically-integrated species concentrations within each of the M pixels, 
indicated by m(k), at each time, t^, from an altitude of h n> v(t^,x^ ,x^ , 

h n ) represents measurement errors. 


Some comments concerning the measurement configuration shown in Fig. 1 
are in order. Ordinarily remote sensing from an airborne platform would be 
carried out at a single altitude. In such a case, it is not possible to esti- 
mate the concentration distribution between the platform and the ground based 
only on the integral of the concentration. Sakawa [8] and Koda and Seinfeld 
[9] have shown that in problems of this nature it is impossible to estimate 
the state uniquely based on integrated measurements from only a single sensor 
position since the required distributed parameter observabil ity condition does 
not hold. Therefore, the estimation of species concentration distributions 
necessitates traverses over the region at different altitudes. From a practi- 
cal point of view this requirement restricts this type of monitoring to air- 
craft platforms, which, for purposes of measuring air pollution, are the most 
useful. Considering that atmospheric concentration distributions change gradu- 
ally and that airplane speeds are fast, the configuration sketched in Fig. 1 
implies that repeated measurements at several altitudes are possible using only 
one airborne platform. 


In order to represent (7) more compactly we introduce the following 


notation: 


x(m(k)) ■ (x* (k >. x™ (k) ) 


m(k) 




m(k)' 3 


(x,) f x~ < h 


3 • 1 _ 

- n 


*1 y H n 

3 n 


J (t k ,x 3 ) 


9l(k) (x 3 5 * * ‘ ? 


5 M(k) (x 3 } 




J( V x 3' 


u(t k , x(l ( k ) ) , x 3 ) I 
u(t k , x(2(k)), x 3 ) 


|^u(t k , x ( M ( k ) ) , x 3 ) 


^ J ( V X 3^ 


L J ^V x 3^ 


f" Z l(k) ( V n) l 


Z(t k ,n) = 




Z M(k)^k’ n ^ 


Z(t L 


Z(t k ,l) 




Z(t k ,N) 


v(t k ,n) = 


v(t k , x(l (k ) ) , h n ) 


v(t k , x(M(k ) ) , h n ) 
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v(t k) 


v(t k J) 


v(t k ,N) 


Then (7) can be represented compactly as 


Z <V * / 0(t k .* 3 )u (* 3 )d* 3 ♦ v{t k ) 


( 8 ) 


We assume that v(t k ) is independent of w(t,x) and u Q (x) and is a white Gaussian 
process with statistics, E [ v ( t. ) ] = 0 and E[v(t k )v'(t e )3 1 R(t k )5 kf , where 
denotes the transpose operator and R(t k ) is an MNxMN positive-definite matrix. 

The problem considered here is to estimate u(t,x) over D on the basis of 

the measurement Z(t ), o * 0, 1 k. The novel aspect of this problem from 

t.he point of view of distributed parameter estimation arises because of the 
scanning and vertically integrated nature of the measurements. In what fol- 
lows, we use k instead of t^ as long as there is no ambiguity. 
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III. Estimation Problems and Wiener-Hopf Theory 

Let us denote the estimate of u(t T ,x) based on the observation data 
Z(t ), o = 0,1,..., k by u(t T ,x/t k ) which is given by the following linear 
transformation of Z(tJ, o = 0,1,..., k, 

k 

u(t,,x/t k ) = Y F(t ,x,t )Z(t ) (9) 

o=0 

where F(t,,x,tJ is an unknown MN-dimensional row vector called the estimation 
kernel function. When there is no ambiguity, we write (9) compactly as 

k 

u(T,x/k) = Y F(t,x,t )Z(t ). (10) 

oSfs T 0 0 

Furthermore, we denote the estimation error and error covariance functions by 
u(t ,x/t k ) and P (t_,x,y/t k ) , respectively, where u(t.,x/t k ) = u ( t T , x ) - u(t T ,x/t k ) 
and p (t_,x,y/t k ) = F. [ D ( t^ x/ t k ) u ( t T ,y/ t k ) ] . The estimate u(t T ,x/t k ) that mini- 
mizes J(u) = E[u(t_,x/t k ) ] is said to be optimal. Note that by using 

A 

P(t.,x,y/t k ), J(u) can be rewritten as J(u) = P(t T ,x,x/t k ) . 

To clarify the differences between the prediction, filtering, and smooth- 
ing problems, we express F(t ,x,t„) differently for each problem as follows: 

(i) Prediction (t > t k ) 

k 

u(t,x/t k ) = Y A(t,x,t )Z(t ) . 
c=0 


(ID 
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(ii) Filtering (t T * t k ) 

k 

G(t k ,x/t k) = ^ F (t k »x,t a )z(t o ). d?) 

(iii) Smoothing (t_ < t, ) 

T K 

k 

u(t_,x/t k ) = £ B ( t T , t R , x , t o ) Z ( t^ ) . (13) 

a=0 

Here we use three temporal arguments t , t, and t„ for the smoothing kernel 
B(t_,t, ,x,t_) since these parameters should be changed according to the measure- 
ment data acquisition time. Then the following theorem can be proved similarly 
to that of [6] for the continuous-time observation case. 

[Theorem 1] (Wiener-Hopf Theorem) 

A necessary and sufficient condition for the estimate u(t, ,x/t k ) to be optimal 
is that the following Wiener-Hopf equation holds for = 0,1,..., k and 
x t D = D _ 3D, 

k 

H F(t,x,t )E[Z(t )Z'(tJ] = E [u(t ,x)Z'( t ) ] , (14) 

o=0 1 *■ c 1 ; 

or equivalently, for £ = 0,1,..., k and x € D, 

E[0(t T ,x/t k )Z'(t_)] = 0. (15) 

[Corollary 1] (Orthogonal projection lemma) 

The orthogonal i ty condition, E[u(t_ ,x/t k )u( t^ ,y/t k ) ] = 0, x,y € D, holds where 

t is any time instant such as t < t. , t = t, , or t > t. . 
n J n k n k n k 


- 10 - 

[Proof] Multiplying each side of (15) by F'U^.y.tJ and summing from 
; = 0 to ; = k yields 

k 

E[0(t.x/t. ) £ Z-(t_)F-(t ,y,t)) = 0. 

T K ~ = q ^ MS 

Using (9) in the above equation yields the desired relation completing the 
proof of the corollary. Q.E.D. 

[Lemma 1] (Uniqueness of the optimal kernel) 

Let F ( t ,x,t ) be the optimal kernel function satisfying the Wiener-Hopf 

T O 

equation ( 14) and let F(t_,x,t_) + F,(t ,x,t ) be also the optimal kernel func- 

i C L TO 

tion satisfying the Wiener-Hopf equation (14). Then it follows that F A (t^,x,t c 
o = 0,1,..., k and x{ D, i.e. the optimal kernel function is unique. 

In order to consider the prediction, filtering, and smoothing problems, 
separately, we rewrite (14) using the notation of (li) - (13). 

(Corollary 2 ) The Wiener-Hopf equation (14) is rewritten for the prediction, 
filtering, and smoothing problems as follows: 

( i ) Prediction ( t > tj, ) 

k 

H A(t,x,tjE[Z(t c )Z'(tJ] = E[u(t,x)Z'(t )] (16) 

o=0 ^ 

for c = 0,1,..., k and x € D. 

(ii) Filtering (t_ = t^) 

k 



ORIGINAL PAOi 
OF POOR QUALITY 
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-li- 


fer C ■ 0,1 , k and x *• 5. 

(ill) Smoothing (t, < t^) 
k 

£ BU..t k .x,t,)E|2(t,)2- 
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(Ul • Elu(t.x)2'(0) 


for C ■ 0,1 


k and x • f>. 
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IV. Derivation of the Optimal Prediction Estimator 

In this section we derive the optimal prediction estimator by using the 
Wiener-Hopf theory in the previous section. 

[Theorem 2] The optimal prediction estimator is given by 

3u(t,x/t.) 

It = L x u(t ’ x/t k } ’ t>l k (19) 

r c G(t,C/t k ) = S(t,0, ^ € 3D. (20) 

[Proof] Differentiating (16) with respect to t and substituting (3) yields 

3A(t,x,t ) 

L rz — E[Z(tJZ'(tJ] = L E[u(t,x)Z'(t )] 

c=0 ?t c x c 

where the independence of w(t,x) and Z(tJ is used. Substituting (16) into 
the above equation yields 

k 

H F.(t,x,t)E[Z(t )Z'(tJ] = 0 

TTn ^ 0 C i. 


where 


3A(t,x,t ) 

F A (t,x,t o ) = ^ L x A(t,x,t o ). 


From Lemma 1 we have 


3A(t,x,t ) 

It " = L x A(t,x, V‘ ^ 


I 


-13- 


ORIGINAL PAGE IS 
OF POOR QUALITY 


* 


Differentiating (11) with respect to t and substituting (21) yields (19). 
Since the forms of 1^ and S(t,S) are known, the predicted estimate u(t,x/t k ) 
also satisfies the same boundary condition (4). Q.E.D. 

[Theorem 3] The optimal prediction error covariance function P(t,x,y/t k ) is 
governed by 


3P(t,x,y/t.) # v 

— — = (L x + L )P(t,x,y/t k ) + Q(t,x,y) , (22) 

r r P(t,C,y/t k ) = 0, K * 3D. (23) 

[Proof] From (3), and (19) we have 

3u(t,x/t.) 

— = L x Q(t,x/t k ) + w(t,x) (24) 

and froT (4), and (20; 

: r u(t,0/t k ) = 0, Ct 3D. (25) 

Differentiating the definition of P with respect to t and using (24) yields 

9. (t.x.y/t^ t L y )P(t,x,y/t k ) + Kt.x.y) 


where 

Z(t,x,y) = E[w(t,x)u(t,y/t k )] + E[u(t,x/t k )w(t,y)]. 
Let the fundamental solution of L x be G(t,o,x,y), where 
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= L x 6(t.o,x,y). 

r r 6(t,o,C,y) = S( t.o , 5 € 3D 

G(c,c,x,y) = 6 ( x-y ) . 


Then u(t,x/t k ) of (24) can be represented in terms of G(t,o,x,y) as follows, 


u(t,x/t k ) = / G(t,t k ,x,a) u(t.,a/t k )da + n G( t,c,x,a)w(r,a)dado. 


\ D 


(26) 


Substituting (26) into I(t,x,y) and using (6) yields I(t,x,y) = Q(t,x,y). 
Multiplying each side of (25) by u(t,y/t k ) and taking the expectation yields 
(23). Q.E.D. 

[Corollary 3] The optimal prediction estimate u(t,x/t k ) and prediction error 
covariance function P(t,x,y/t) can be represented as 

u(t,x/t R ) = / Glt,t k ,x,a)u(t k ,a/t k )da (27) 


and 

P(t,x,y/t.) = // G(t,t, ,x,a)P(t, ,a,S/t, )G(t,t. ,y,B) dad£ 

k D d K 

+ j JJ G(t,o,x,a)Q(o,a,£)G(t,c,y,6) dadBda. (28) 

t. D D 
k 

[Proof] It is clear that (19) and (22) possess unique solutions. Differen- 
tiating (27) and (28) with respect to t yields (19) and (22), 
respectively. Since (19) and (22) have unique solutions, (27) and (28) are 
those solutions. Q.E.D. 
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V. Derivation of the Optimal Filter 

In order to derive the optimal filter by using the Wiener-Hopf theorem 
for the filtering problem, we represent the solution of (3) in terms of the 
fundamental solution G(t, o,x,y) as 


u(t. 


»*) - f G( t, , -I ,t. ,x , 


a)u(t, ,a) da 


t k+1 r 

J J G(t k+ 1 ,n,x,a)w(n,a) 


dadn 


\ D 


(29) 


and 


u t k+1 ( * 3) = l Wl’W^^V 0 ^ d: 


r k+1 r 

+ J J GJ t, + , ,n»x 3 a)w(n,a) dotdn 


l k D 


(30) 


where 


t k+ ^ , r ,X2»a) = 


rlt n l(k + l) l(k+l) „ „» 


G( t k+ ^ » r ' ’ X- | 


M( k+1 ) M( k+1 ) 


,*2 ’ x 3 * 


(31) 


Frorr (17) we have 


F(t k+l’ x ’W E[Z(t k+l )Z ' ( V ] 


+ X F(t k+ 1 ,x,t o )ElZ(t c )Z'(tJ) = E[u(t k+ 1 ,x)Z'(t ; )l 


(32) 


for : = 0,1,. , k+1. 
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From (29) and the independence of Z(t^), c = 0,1,..., k and w(n,x), 

t. < n < t. it follows that 
k — k+1 

E[u(t k+1 ,x)Z - (t ; )] = / G(t k+1 ,t k ,x,a)E[u(t k ,a)Z'(t ; )] da. 


Using the Wiener-Hopf equation (17), we have 


E[u(t k+1 ,x)Z'(t ; )l = / G(t k+1 ,t k ,x,a) II F(t k ,a,t 0 )E[Z(t 0 )Z-(t c )]. 


On 


the other hand, from (8) and the whiteness of v(t k+ j), we have, for 


- V 


f ^ 

E[Z(t k+1 )Z'(tJ] = J J(t k+1 ,x 3 )E[u (x 3 )Z'(t )] dx 3 . 

o k+1 

Substituting (30) into the above equation and using the independence of 
Z(tJ, t. < t k und w(n,a), t k < n < t k+1 yields 

E[Z(t k+1 )Z'(tJ] = / J(t k+1 ,x 3 ) / G M (t k+1 ,t k ,x 3 ,a)E[u(t k ,a)Z > (t ; )] dadx. 


Again, we use the Wiener-Hopf equation (17) in the above equation and 


E[Z(t )Z'(t )] dadx, 

J 


h k 

E[Z(t k+ i)Z (t ) ] = / J(t k+1 ,x 3 ) / G M ft k+1 ,t k ,x 3 ,a) I^F(t k ,a,t c ) 

n 0 = U 


( 34 ) 


Substituting (33) and (34) in (32) yields 
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o=0 


^V X ’V E z <v r <v 


where 


^ A (t k » x .t a ) = F(t k+1> x ,t k+ i) f J ( t | c +j » x 3^ / G M^ t k+r t k’ X 3’ 0l ) F ^ t k ,a ’ t '^ dadx : 


+ F ^ Vr x, V ’ ^ G ^ t k+r t k ,x ,a)F(t k ,a,tJ da . 


k o 


Then from Lemma 1 we have F (t. ,x,t ) = 0, and we have the following lemma. 

u K C 

[Lemma 2] The optimal kernel function F(t k+ j,x,t„) of the filter is given by 


F ( t k+1 ,x,t o) = | G( tfc+l *t k .x,o.)F( t k ,a,t c ) 




r 

- F (t k+1 ,x,t k+1 ) J t | C+ i - x 3) J [ dadx 3* 

0 D 


[Theorem 4] The optimal filtering estimate u ( ^ ’ X/ ^ 1 ^ g ’ ven by 

u 'tk+r x/t k+l) = u ^k+l ,x ^k^ + F ^k+1 ,x **k+l^ v ^k 4 '» ^6) 


v(t k+ i) E Z(t k+1 ) - / 'J( t k+r x 3 ) u t k+1 (x 3 /t k ) 
u(t c ,x/t 0 ) = u Q (x), 


-U ( t k+ j , c/ t k+ j ) - r <E 9D 


dx. 


(37) 

(38) 

(39) 


where 
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U* ( X-/ ) 

L k +1 J * 


“( t k+r x l^ k * 1 ^ x 2^ k+1 ^ ,x 3 /t k^ 

u(t x M(k+1 > x M(k+1) X /t ) 

ul Vl’l ,x 2 * x 3 /l k'_ 


[Proof] Using (12) and (35) yields 


u(t k+r x/t k+l) = F ^ t k+r x,t k+P Z ^ t k+l^ 


* 

1 G(t k+r t k .x,a) YL Fd^.a.tJZUJ dot 


k*“* V 


F ' t k+r x,t k+l) / J(t k+r x 3 ) I G M^ t k+l ’ t k ,x 3 ,0t ^ ^ F ^ t k* a,t o^ Z ^ t o^ dadx : 


Then frorr. (12) and (27) we have 


' 

U ^k+1 ’X/^k+l) - ^ ( ^k+ 1 * ^ k * x * ^) u ( ^ k * ^ ^k ) ^ 


,h r 

F ( t k -*-i ,x ’ t k+i ) ( z ( t k+i ) " i J(t k+r x 3^ l Si^k+rV^’^^V^ V dadx : 


u(t k+r x/t ! ) + F ^ t k+r X,t k+l- v ^ t k+l^ ■ 


Since the initial and boundary conditions are clear, the proof of the theorem 
is complete. O.E.D. 

To determine the optimal kernel function F(t k+1 ,x,t k+ ^) , we introduce 
the following notation, 

P M (t_,x,y 5 /t k ) = [P(t,,x,y 1 ' k ^/t k ),..., P(t,,x,y M ^ k Vt k ) ] (40) 


and 
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WWW 


P(t T ,xl (k) ,y 3 /t k ) 


P(t T ,x M ^ k ^,y 3 /t k ) 

P(t T ,x 1 (k) ,y 1 (k) /t k ) P(t T ,x 1 (k) ,y M(k) /t k ) 

• • 

P<t T .»" (; V ,k >/t k ) P(t T ,x M ' k >. y H < k )/t k )_ 


(41) 


. m(k) , m(k) m(k) x m( k ) _ / m( k ) w m ( k ) v \ m - ? m 

where x v - (x^ x^ , x 3 ) and y - (yj , y 3 » yjl* m ~ !»*»•••» n . 

From the definitions of P^( t_ , x ,y 3 / 1^ ) and Pf^^ * x 3*- y 3^ t k^ ^ ^°^ ows that 


P M (t V X,) W = E[u(t T ,x/t k )u t (y 3 /t k )l 


(42) 


and 


P Mfr t T * x 3 * y 3 / t k ^ ~ E ^t ( x 3 /t: k^t ^3^*k^ 


(43) 


where 


u t (x 3 /t k ) * u (X 3 ) - u t (x 3 /t k ) 


(44) 


and 


u t •W = 


u(t. ,x|' k >/t k ) 

• 

M( k) . X 
u ( 9 x / 1 ^ ) 


(45) 


Furthermore, we define the covariance matrix of the innovation process v(t k+ ^) 


by r(t k+ 1 /t k ) 

= E[v(t k+1 

r <WV 

,h ,h 

-f f J(t 

0 0 


* R( W- 


( 46 ) 
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[Theorem 5] The optimal filtering gain function F( ,x , is 
given by 

F( Vr x, W = / p M (t k*r x,y 3 /t k )J 'Vr y 3^ dy 3 r ^ t k+: /t k ^ 

u 

[Proof] From the Wiener-Hopf equation (17) we have 

F(t k+ 1 ,x,t k+ i)E[ z ( t k + i)Z (t k+1 )l 


r\ 

+ T, F!t k 4 l .x,t_)E|Z(t c )Z'(t k+1 )] ■ EIu(t k+ 1 .x)Z'(t k+ ,)l 


0 = 0 


k+1 ' v k+1' 


Substituting (35) into the above equation yields 


Flt k4r x ’Vi )E i ( z (t k+r " / J(t k+r x 3 )u t. +1 (x 3 /t k )dx 3 )z (t k+i' ] 

o k+1 


E [ ( u ( t k+1 . x ) - u(t k+1 ,x/t k ))Z'(t k+1 )] 


Using (8) ana the orthogonality condition of Corollary 1 yields 


Elull .x/t )Z'(t )] = / E!u(t ,,x/t k )u' (x 3 )]J'(t ktl .x 3 ) dx 

o k+1 


r 

J V t k+r x,x 3 /t k' J (t k+r x 3 ) dx : 


and 


f h ,h 

E ( v (t k4 i)Z'( t k+ 1 ) ] = J J J l t k+l’ x 3^ P MM (t k+r X 3 ,y 3 /t k^ J ^k+l ,y 3^ dx 3 dy ) 

o o 


n*»rr 
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r h 

^Vl^’W^WV = i P M lt k+r x,x 3 /t: k^ (t k+l ,x 3 } dx 3 


and the proof of the theorem is complete. Q.E.D. 

[Theorem 6) The ootimal filtering error covariance function Pvt k+ r x .y/t k+1 ) 
is given by 



p(t k+r x 


P M ( Vl’ x ’ x 


•y / Vi ) = p(t k+r x,y/t k ) 

3 /l k^ J ^V+l’V ^ t k*l / V‘^ t k*l 


* y 3^ P M^ ^k-H * y ,y 3 /t; k ^ dx 3 dy 3 


P( V x ’ y/ V = p 0 (x ’ v) ’ 

: r P(t k+1 ,C. y /t k ) =0, l, «* 3D. 

[Proof] From (3) and (36) we have 

u(t k^l ,x/t k+l^ = ^k+r^V " F( Vr X,t k+l )v( Vl 

and from (4) and (39), 

V^kT^W = °' 3D ' 

Using the independence of v(t k+1 ) and u(t k+ ^ .x/tj.) or G(t k+ j 

p (t k+ i ,x, y /t k+ i) = E[u(t k+ j,x/t k+ j)u(t^ <f |, y /t k<f i) ] 


(49) 

(50) 

(51) 

) (52) 

( 53 ) 

,y/t k ) yields 
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p (t k+1 ,x,y/t k ) + p ( t k+ i’ x » t k+1 ) E fv(t k+1 )v (t k+1 )]F 

( t k+l ,y,t k+l ) 


F ( t k +l’ x,t k + 1 ) i J ^ t k+l fX 3^ Iu t. ./V V u ^k+r y/ V J dx 3 
c k+ 1 

r h 

] E[u( t k+ i, x /t )( )D t (y 3 / t k )] J (t k+l’ y 3 )dy 3 F ^ t k+l* y • t k+l ^ ’ 


Using (40) and (47) yields 


P(t k+1 ,A,y/t k+1 ) = P(t k+1 ,x,y/t k ) 


/ / P M' t k+l ,x,x 3 /t k )J (t k+l’ x 3 ): 1 ( t k+l /t k ^ 


J, 'Vr y 3' dx 3 dy 3* 

A — 

Since tne initial value u(t Q ,x/t 0 ) is equal to u Q (x), it is clear that 
P(t o ,x,y/t Q ) = P (x,y). Multiplying each side of (53) by u( t k+ j ,y/t k+1 ) and 
taking the expectation yields (51). Q.E.D. 
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VI. Derivation of the Optimal Smoothing Estimator 

In this section we derive the optimal smoothing estimator by using the 
Wiener-Hopf theory. 

[Lemma 3] The optimal kernel function B(t T ,t^ + ^x,t_) of the smoot^ng estima- 
tor is given by 


B ( t T »t k+1 ,x,t c ) = B(t T ,t k ,x,tJ 


® ' ^k+1 ’* ’ ^k + l ^ J J(^k+1 .* 3 ) f »t k . X 3 ,a) F ( t k ,a , t^ ) dadx^. (54) 


[Proof] From the Wiener-Hopf equation (18) for the smoothing problem we have 
k+1 

X B(t ,t k+1 ,x,tJE[2(t )Z'(tJ] = E[u(t,x)Z'(tJ], (55) 

r=u 1 x i c s t ? 

; = 0,1 k+i 


and 

k 

^ B(t_ ,t ,x,t )E[Z(t )Z'(t ) ] = E[u(t .x)Z'(t)], (56) 

k ~ a q 1 c, 

0- u 

; = o,i...., k. 

Subtracting (56) from (55) yields 

B ^ t T ,t k+r x ’ t k+l^^' t: k+l^ Z ^5^ 
k 

+ Swvvi-.v ■ B (vv x,t m ))E[Z(t a )r( v ] = °* 


From (8) and (17) wc have 
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EIZ(t k+ ,)Z'(t )] = / J(t k+1 -x,)/ G M {t |c+ 1 ,t k ,x 3 c I )Etu(t k ,a)Z'(t c )] dadXj 
^ o 0 

h - k 

= J J ' t k+l’ x 3^ ( ’M^ t k+l* t k’ X 3’ a ^ ]^ F (V a »V E[Z ^V Z dadx 3‘ 


Then it follows that 


1 Q ? A-V x * t o )EIZ(t o ,Z ' ( V 1 


= 0 


where 


^VV'-V * B(t T ,t k+r x,t 0 ) - B(t T ,t k ,x,t 0 ) 


* B l t -’ t k+l’ x,t k+l ) / J!t k*l’ x 3' / G H (t k+r t k’ x 3’ a * F(t k’ a,t o* 
‘ 0 u 


Since it is clear that B(t T ,t k ,x,t,) + F,(t T ,t kl x,t 0 ) also satisfies the Wiener- 
Hopf equation (18), from Lemma 1 F , ( t^ , t k , x , t^ ) = 0, o = 0,1,..., k. Thus, 
the proof of the lemma is complete. Q.E.D. 

[Theorem 7] The optimal smoothing estimate u(t .x/t k+1 ) is given by 

u(t T ,x/t k+ ^) = u(t T ,x/t k ) + ^^ T »^k+i ,x, ^|(+i^ v ^k+l^ 

r^u(t T ,c/t k+1 ) = S(t,0» l 3D - ( 58 ) 

[Proof] From (13) it follows that 

u ( t_ , x/ ) = B ' t T » t ( c +i’ x ’ t k+l^ Z ^ t k+l^ 
k 

+ X ^Wl’^O^ V* 

c=0 
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Substituting (54) into the above equation yields 


u ( t T » x / t k+ i ) = B( VVr x,t k+i )v(1 W 


£ B(t T ,t k ,x,t 0 )Z(t a ) 
o=0 


and substituting (13) into the above equation yields (57). Since we have no 
additional information about the boundary value of u(t T ,x) except for S(t,,5). 
we have (58). Thus, the proof of the theorem is complete. Q.E.D. 

[Theorer 8] The optimal smoothing gain function B(t T ,t k+ - ,x,t^ + j) is given by 

_ i 

B (t_,t k+ r x ,t k+ i) = J N(t.,x,x 3 /t k+1 )J'(t k+1 ,x 3 )dx 3 r (t k+1 /t k ) (59) 

0 

where 


N(t T ,x,x 3 ./t k+1 ) = / M(t T ,x,y/t k )G'(t k+1 ,t k ,x 3> y) dy 


and 


Mvt..x,y/t k ) = E[u(t_,x/t k )u(t k ,y/t k )]. 
[Proof] From the Wiener-Hopf equation (18) we have 


+ 


B ( tj , t k+ j » x »t k+ j)E[Z( t k+ ^ ) Z ( ^ k+ j ) ) 
k 

£ B(t T ,t k+1 ,x,t o )E[Z(t c )Z'(t k+1 )] = E[u(t T ,x)Z'(t k+1 )]. 

o=0 


(60) 


(61) 
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Substituting (54) into the above equation yields 


B (VVr x,t k+l )E[vl Vl )Z,(t k + l )1 = E[u( V x/t k )Z (t k+l }1 

On the other hand, froir (27) and (29) 

- / G(t kM ,t k .x,y)u(t k ,y/t k ) dy 

A+if 

+ I J G(t. + 1 ,n,x,y)w(n,y) dydr,. 

*k D 


( 62 ) 


Then we have 


E[Q(t .x't. )Z'(t ,) J * / / MU..x,y/t k )G-(t k+r t k ,x 3 .y) dy J'(t k4 ,.* 3 ) d » 3 

K 0 D 

= / N(t. ,x,x 3 /t k+ 1 )J'(t k+ 1 ,x 3 ) dx 3 . 


Substituting (48) and the above equation into (62) yields (59). Thus, the 

proof of the theorem is complete. Q.E.D. 

Let us now derive the equation for M(t T .x.y/t k+1 ). Using the orthogonality 

condition of Corollary 1 yields 

Mit T .x,y/t k+1 ) * E[u(t T .x)a(t k+1 .y/t k+1 )1. (63) 

Substituting (52) into the above equation yields 


Mv 


t..x.y/t k+ .) - / Sit k+ 1 .t k .y,a)M(t T ,x,a/t k ) da 

Hi q 


- / N(t ,x,x 3 /t k+ 1 )J'(t k+ 1 ,x 3 ) dx 3 F'(t k+r y,t k+1 ). 
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From (4) and (58) we have 


r 5 u(t T , 5 /t k+1 ) - o, 


E, T. 9D. 


(64) 


Multiplying each side of the above equation by d(t k+ j,y/t k+1 ) and taking the 
expectation yields r £ M{ t T . C ,y/t k+1 ) ■ 0. « € 30. Thus, the following theorem 
holds. 

[Theorerr 9] M(t ,x f y/t R+1 ) is given by 


M(t_.x,y/t k+ j) = / G(t k+1 ,t k ,y^)M(t T ,x,a/t k ) 


da 


(65) 

( 66 ) 
(67) 


- / N(t T ,x,x 3 /t k+ 1 )J'(t k+ 1 ,x 3 ) dx 3 F '( t k+ r y,t k +l) ’ 
u 

M( t_ ,x ,y/t_ ) = P(t T ,x,y/t_) , 

r E M(t T .;,y/t k+1 ) * o, fr3D. 

It remains to derive the equation for the optimal smoothing error covariance 
function P(t_ ,x,y/t k+1 ) . From (57) we have 

D(t.,*/t k+1 ) • D(t,,x/t k ) - B(t T ,t k+1 ,x,t k+ i)v(t k+I ). (68) 

[Theorem 10] The optimal smoothing error covariance function P(t T ,x,y/t k+ , ) 
is given by 

P(t,x,y/t. +1 ) = P(t ,x,y/t k ) 


- f j N(t .x,x 3 /t k+1 )J-(t k+1 .x 3 )r A ( J( t k+l ,y 3 ^ 

o o 


- 1 , 


N(t ,y f y 3 /t k+1 ) dx 3 dy 3 , 


(69) 
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r.P(t T ,f,y/t k+1 ) = 0, i € 3D. (70) 

[Proof] Frorr ( 68 ) we have 

P(t r ,x,y/t k+1 ) = E[u(t T ,x/t k+ 1 )u(t T ,y/t k+1 )] 

= P(t T ,x,y/t k ) + B (t T .t k+ 1 ,x,t k+ 1 ), (t k+ 1 /t k )B (^ T .t k+ 1 ,y,t k+1 ) 

‘ B ^ t x’ t k+r x ’ t k+l^^ v ' 

- E[Q(t T ,x/t k )v'(t k+ 1 )]B-(t T ,t k+ 1 ,y,t k+1 ). (71) 

But we have 

E [ u . t_ ,x/t k )v (t k+ j)] = f J ^ ’ X 3 » a )M( t^ ,x , 0 ./ t k ) 

0 D 

J'(t k+ r x 3 ) d ^ dx 3 
and 

EJv(t k+ 1 )u(t ,x/t k )] = (j J(t k+ 1 ,x 3 )G L (t k+ 1 ,t k ,x 3 ,a)M(t T ,x,a] d^dx 3 . 

0 D 

Substituting the above equations and (47) into (71) yields (69). Multiplying 
each side of (64) by u(t„,y/t k+1 ) and taking the expectation yields (70). Q.E.D. 

[Theorem 11] The optimal smoothing estimator is given by 

k 

u(t ,x/t.) = u(t.x/t.) + 2 B(t ,t ,x,t )v(t ) 

k T r v=T+1 T X X X 


( 72 ) 
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and the optimal smoothing error covariance function P(t T ,x,y/t k ) is given by 
P(t.,x,y/t k ) = P(t T ,x,y/t T ) 


* h h 

- 1 / /'N(t ,x,xj/t,)j'(t t .* 3 )r* 1 (t t /t t . 1 ) 


=: + l 0 0 


J(t ; ,y 3 )N(t T ,y,y 3 /t Jl ) dx 3 dy 3 . 
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VII. Estimation of the Concentration Distribution Downwind of a Continuous, 
Ground-Level Line Source 

There has been much recent interest in the airborne measurement of pollu- 
tant concentrations downwind of sources [9 ] - [11]. Here we wish to consider 
a hypothetical, but realistic, situation in which an aircraft with a downward- 
looking instrument, such as for example the JPL Laser Absorption Spectrometer 
[12], is flown at different altitudes downwind of the source, and total species 
burdens are measured at a series of downwind distances. 

The steady-state concentration of a species downwind of a continuously 
emitted ground-level line source (e.g. a highway) situated normal to the 
direction of the wind flow is governed by the following form of the atmospheric 
diffusion equation [5] . 


v ly. 
‘1 5x. 


c*X. 


K (x 

w ' 


?u 


'1 l,A 3 ' 
u(0,x 3 ) = u q (x 3 ) 


v VA 3' Bx, I " W(x l’ x 3 ) 


“ U°) TT = x 3 = 0 


?x 


(74) 

(75) 

(76) 


3 


Bx 3 


= 0, 



(77) 


where £ is the constant rate of release. For convenience we will take = 1, 
since vertical variations of this constant are not essential to the estimation 
problem we will consider. If we let t = x^/Vj and x = x 3 , (74)-(77) become 


f>U = £_U + ^ y \ 

at ,2 V ' l ’ x) 

Bx 


u(0,x ) = Ug( x] 


(78) 


(79) 
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t * - ♦«(»), x = 0 

lx ■ °' x * h 


(80) 

(81) 


In this case the measurements Z(t, ) are related to the concentration 

In 

u(t k ,x) by (8), 


f h 

Z(t k ) = J J n (t k ,x)u(t k ,x) dx + v (t k ) (82) 


where the instrument kernel function will be taken to have the form, 

r 


J (t k ,x) = 


x < h 
— n 


x > h 


n = 1,2,..., N 


(83) 


The theory developed in the prior sections can be applied directly to this 
problem, and the optimal filter and smoother are given in Table 1. The pre- 
diction, filtering and smoothing algorithms were applied to hypothetical data 
generated by solving ( 74 ) - ( 77 ) and forming Z(t^) from (82), using noise 
processes w(t,x) and v(t k ) with prescribed properties. The algorithms were 
applied to estimate the concentration distribution u ( t k , x ) as a function of 
height x at several downwind distances, t^ , t^,... based on measurements taken 
at one to four elevations. It is of interest to study the behavior of the 
estimates as a function of downwind distance and of the number of elevations 
at which data are simultaneously taken. Values of a’l parameters used in the 
calculation are given in Table 2. 
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Figs. 2-4 show selected results of the application of the filtering and 
smoothing algorithms to the synthetic data of this example. Fig. 2 shows a 
comparison of the true concentration distribution u(tj,x) and the filter esti- 
mates, u(tj,x/tj) based on two and four measurement elevations (t^ = 0.0002). 

As expected, the profile estimated on the basis of four measurement elevations 
is superior to that based only on two altitudes. Fig. 3 shows similar results 
at tg = 0.0082. The filter estimate based on n = 4 virtually coincides with 
the actual concentration distribution. The performance of the smoothing algo- 
rithm is illustrated in Fig. 4, in which the true concentration u(t T ,x) is com- 
pared with the filter estimate, u(t_,x/t„), and the smoothed estimates, 
u(t ,x/to), and u(t,,x/t^), with t T = 0.0002, tp = 0.0012, and t^ = 0.0032. 

Table 3 gives the trace of the filtering error covariance matrix, P(t,x,x/t), 
for the four measurement configurations at three downwind distances t. As expec- 
ted, the trace decreases as the number of measurement elevations is increased 
from 1 to 4. 

VIII. Conclusions 

Filtering and smoothing algorithms for the processing of remote sensing 
data on atmospheric species concentrations have been derived using Wiener-Hopf 
theory. The algorithms were applied successfully to estimate concentration 
distributions from a hypothetical ground-level line source of material (e.g. a 
highway) based on remote sensing data taken from several elevations at a number 
of points downwind from the source. Although there has been increasing interest 
in the remote sensing of airborne concentrations, a data set sufficient for ap- 
plication of the theory developed in this paper does not yet appear to exist. 
Nevertheless, it is hoped that the availability of the algorithms developed here 
will facilitate processing of remote sensing data in conjunction with mathematical 
models of air pollutant behavior. 
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Table 1. Optimal filter 


Filter 


and Smoother for Line Source Application (Continued) 


Smoother 


ii) (Continued) 


- f J(t T ,f»)N(t T ,x,a/t k )daF'(t k ,y,t k ) 

r(l. ../t. ) = / h / hj( S+l •*> p< ‘ktl ’* ,y/t k' J( t k*r v) ‘ 


dxdy + R(t k+ j) 

P(t. Xl ,x,y/t. ..) = P(t k+r x,y/t k ) 


M(t i ,x,y/t i ) = P(t T ,x,y/t 7 ) 


/Y J|t W“ ,P(, W'''" /l k ) 


r ' 1 <W t k )P(t k*r B,y/t k ) d ' ,,)P 


s 


ORIGINAL PAGE IS 
OF POOR QUALITY 


origin* 1 - 

0F POOR QUALITY 


Table 2. Parameter Values Used in Line Source Estimation 

Example 


Truncation number N = 5 


Measurement time t k+1 = t k + 0.001, k = 1,2,3,-’* 
where t-j = 0.0002 

Fixed- point time for smoothing t = 0.0002 
Constant rate of release <f = 0.3 

Measurement points h n = n/4, n = 1, 3,4 

Initial values and noise covariances 

E[Uq(x)] = Uq(x) = ; U°0.(x) 

N . 

Cov[u Q (x),u 0 (y)] = P Q ( x,y ) = J p° ^ (x)C. (y ) 


Cov[w;t,x),w(s,y)] = Q(t,x,y)f(t-s), Cov[ v( t k ) , v(t n ) ] = R ( t k )- kn 
N 

Q(t,x,y) = l q ii 4'.(x);.(y), R(t fc ) = diag[r ] ,r 2 ,r 3 ,r 4 ] 
i=l 


C i (x) 


1 

>7 c o s ( i - 1 )-x 


i = 1 

i > 2 
i > 1 


I 

i i 

1 

— 

2 

3 

4 

— 

5 

hr 

3.0 

1.0 

0.03 

0.003 

0.0003 

pr 


0.1 2 

0.01 2 

o 

§ 

no 

0.0001 2 

*11 

i 

i 

0. 5 L 

C.25 2 

0.125 2 

0.0625 2 

1 r. 

i i 



0.1 2 

0.07 2 

0.05 2 

0.03 2 





- a 





****■ 1 1 ii n 



Table 3. 


Trace of the Filteri 


Measurements 

t = 0.0002 

4 point 


(h-j » » h^ » h^ ) 

0.2405x10' 

3 point 


(h^ * h£ » h^ ) 

0.3441x10' 

2 point 


(h ] ,hj 

0.1678 

1 point 


(h ] ) 

0.6700 
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Figure Captions 


Fig. 1. Remote Sensing Measurement Configuration Considered in This 
Work. 

Fig. 2. Comparison of True Concentration u(t-|,x) and the Filter Esti- 
mates u(t^,x/t^) based on 2 and 4 Measurement Elevations. 
t ] = 0.0002. 

Fig. 3. Comparison of True Concentration u ( tg , x ) and the Filter 

Estimates u(tg,x/tg) based on 2 and 4 Measurement Elevations. 

tg = 0 . 0082 . 

Fig. 4. Comparison of True Concentration u(t ,x), the Filter Estimate 
u(t. ,x/t T ) and the Fixed Point Smoothing Estimates u(t T »x/t 2 ) 
and u(t T ,x/t 4 ). t, = 0.0002, t 2 = 0.0012, t 4 = 0.0032. 
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SUMMARY AND CONCLUSIONS 

The object of this research grant was to initiate an evaluation of the 
analysis of remote sensing data on pollutant concentrations in the troposphere. 
Remote sensing measurements of pollutant concentrations are becoming increas- 
ingly important in understanding the transport and transformation of pollutants 
over moderate to long distances in the atmosphere. Traditionally such data 
have not been analyzed beyond the point of constructing mass fluxes and total 
budgets over a region. The question studied in this research grant w s that 
of the further analysis of such data, particularly when one has a mathematical 
model available. The specific problem then is to see how typical remote sens- 
ing data can be used in conjunction with a mathematical model to extract addi- 
tional information about the pollutant behavior in the region being studied. 

The essential problem is one of estimation, that is, of using the typical 
remote sensing data to determine full concentration distributions. Once full 
concentration distributions are available, one can then assess the mechanisms 
of the process through the mathematical model. The first step in the research 
was to look theoretically at the question of the minimum amount of data needed 
to reconstruct a concentration distribution from finite data typical of those 
collected in remote sensing. Chapter I of this report presents a development 
and derivation of a condition of reconstructability, namely rigorous conditions 
that can be applied to a data sampling program to determine whether it will be 
possible to estimate a species concentration distribution from such measure- 
ments. Chapters II and III of this report are then devoted to the development 
of a numerical algorithm that will process the data to produce concentration 
distrioulion estimates in the cases when the data are a priori reconstructable. 

Perhaps the most important result of this study is the indication of the 
types of measurement strategies one should follow in remote sensing programs. 


In particular, it appears that the best measurement strategy is to attempt to 
obtain pollutant burdens at a certain location at a number of elevations at 
t mes as close as possible. This strategy is recommended because the vertical 
distribution of pollutant concentrations in the first 1,000 meters of the atmo- 
sphere is a crucial element of a mathematical model of such species. The 
theory and numerical techniques developed in this study will tell one when 
devising a measurement program and monitoring strategy the number of vertical 
levels at which one should make measurements to be able to estimate relatively 
accurately the complete vertical concentration profile of the soecies of inter- 
est. It is anticipated that these results will be of value to those contem- 
plating remote sensing measurement programs of tropospheric species that involve 
measurements at several vertical levels. 


